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FLOWFIELD ANALYSIS FOR SUCCESSIVE 
OBLIQUE SHOCK WAVE-TURBULENT 
BOUNDARY LAYER INTERACTIONS 


By Chen-Chih Sun and Morris E. Childs 
Department of Mechanical Engineering 
University of Washington 


SUMMARY 

A computation procedure is described for predicting the flowfields which 
develop when successive interactions between oblique shock waves and a turbu- 
lent boundary layer occur. Such interactions may occur, for example, in engine 
inlets for supersonic aircraft. Computations have been carried out for axisym- 
metric internal flows at M_^ = 3.82 and 2.82. The effect of boundary layer bleed 
has been considered for the M = 2.82 flow. A control volume analysis is used to 
predict changes in the flow fTeld across the interactions. Two bleed flow models 
have been considered. A turbulent boundary layer program has been used to com- 
pute changes in the boundary layer between the interactions. The results given 
are for flows with two shock wave interactions and for bleed at the second inter- 
action site. In principle the method described may be extended to account for 
additional interactions. The predicted results are compared with measured re- 
sults and are shown to be in good agreement when the bleed flow rate is low (on 
the order of 3% of the boundary layer mass flow), or when there is no bleed. 

As the bleed flow rate is increased, differences between the predicted and mea- 
sured results become larger. Shortcomings of the bleed flow models at higher 
bleed flow rates are discussed. 


INTRODUCTION 

The interaction of an oblique shock wave with a turbulent boundary layer is 
known to induce drastic changes in the boundary layer properties and to cause 
substantial deviation of the supersonic flow field from the predicted inviscid 
flow. This deviation may be of sufficient magnitude to adversely affect the per- 
formance of aerodynamic devices. Suitable methods for predicting the boundary 
layer and the freestream flow characteristics in the presence of such disturbance 
are required by engineers responsible for the design of aerodynamic configurations 
in which shock wave boundary layer interactions occur. 

A control volume method developed by Seebaugh, Paynter and Childs [1], and 
improved upon by Mathews [2], has been used successfully in the prediction of 
the boundary layer characteristics downstream of the interaction with a single 
oblique shock wave. However, in some aerodynamic devices such as mixed compres- 
sion supersonic diffusers, the turbulent boundary layer is subjected to inter- 
actions with more than one shock wave. In this report, a computation procedure 
is described for predicting the flow field which develops when successive inter- 
actions of two oblique shock waves with a turbulent boundary layer occur. In 
principle the method described may be extended to account for additional inter- 
actions. 


Computations have been carried out for axisymmetric internal flows at = 
3.82 and 2.82. Experiments have been conducted at these same Mach numbers 
with the interactions under study occurring at the walls of circular wind tun- 
nels. In the Mach 2.82 study the effect of boundary layer bleed at the second 
interaction site has been considered. The predicted results are compared with 
experimentally observed results. The experimental configurations for which 
the analysis has been carried out are discussed in the section which follows. 
This is followed by a section which gives the details of the computational 
method. 


SYMBOLS 

A = {[(Y-l)/2]Mg2/(TyTg)}^^ 

a = a constant in the wall -wake profile 

B = {(H[(Y-l)/2]Mg2)/(yTg)}-l 

C = a constant in the Law of the Wall (usually equals 5.1) 

= skin friction coefficient, ■^v^/n/ 2 )pgUg 2 
F = entrainment function, see Eq. (B-5) 

(H.|)|^= kinematic shape factor, see Eq. (B-6) 

Ibx = x-momentum of the bleed flow 

K = a constant in the Law of the Wall (usually equals 0.4) 

L = shock wave boundary layer interaction length 
M = Mach number 

mg = boundary layer mass bleed rate 
P = pressure 

R = radial coordinate from tunnel centerline 

Rg = radial coordinate of dividing stream surface separating bleed flow 

from main flow, (see Fig. 3) 

Rg = Reynolds number 

u = velocity in streamwise direction 

■I / 

u = vanDriest's generalized velocity, (Ug/A) arcsin {[(2A2u/Ug) - B]/(b2+4A2) } 
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= friction velocity, (t /p.,) “ 

X = axial coordinate, measured from shock generator tip 

y = coordinate normal to the tunnel wall 

Y = ratio of specific heats 

A = mass flow thickness, see Eq. (B-1) 

Ac = a thickness of freestream flow to allow for boundary layer mass 
entrainment, (see Fig. 3) 

6 = boundary layer thickness 

6* = displacement thickness of the boundary layer 

n = y/6 

9 = momentum thickness of the boundary layer 

V = kinematic viscosity 

n = coefficient of the wake function 
p = mass density 

o = [(Y-l)/2]Mg2/{i+[(Y-i)/2]Mg2} 

T = shear stress 


Subscripts 

e = conditions at the edge of the boundary layer 

w = conditions at the wall 

CO = freestream conditions ahead of the first interaction 


EXPERIMENTAL CONFIGURATION AND INSTRUMENTATION 

The experimental configurations which were used to produce the successive 
shock waves are shown schematically in Figure 1. The Mach 3.82 tunnel had a 
radius of 2.58 cm and a boundary layer thickness ahead of the first interaction 
of 0.43 cm. The shock wave generator was installed on the centerline of the 
tunnel at zero angle of attack. The generator had a 10-degree half-angle coni- 
cal tip which broke to 13 degrees 6.60 cm downstream of the tip. For the Mach 
2.82 tunnel the tunnel radius was 2.60 cm and the boundary layer thickness 
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ahead of the first interaction was 0,40 cm. The shock generator had a 10- 
degree conical tip which broke to 13 degrees 4.94 cm behind the tip. Both 
generators were designed to provide as large a region of freestream flow as 
possible between the first reflected and second incident shock waves while at 
the same time keeping the expansion wave off the downstream corner of the 
second conical surface from interfering with the second interaction. 

About 7.62 cm from the tip of the Mach 2.82 shock generator, or approxi- 
mately at the point where the second incident shock wave reached the wall, two 
rows of thirty-eight 0.132 cm diameter bleed holes were drilled around the 
periphery of the tunnel. The bleed system was operated in a choked flow con- 
dition. With one row of the holes open, the bleed mass flux was about 2.8 
percent of the boundary layer mass flux just ahead of the second interaction. 
With two rows of holes open, the bleed mass flux was 5.0 percent of the boun- 
dary layer mass flux. 

Both tunnels were operated with a steady supply of dry air at 300°K. The 
freestream unit Reynolds number for the M = 3.82 tunnel was 18.5 x 10® per 
meter, that for the M = 2.82 tunnel, 19.0 x 10® per meter. 

Standard instrumentation was used to obtain tunnel wall static pressures 
and boundary layer oitot pressures profiles. Wall static pressures were taken 
at 0.127 cm intervals along the tunnels. Pitot profiles were taken 

in radial increments of 0.0127 cm at eighteen axial stations upstream of, with- 
in, and downstream of the interaction region. Miniature total head tubes, 
flattened to a dimension of 0.024 cm high by 0.066 cm wide were used for the 
pitot profiles. Velocity profiles upstream of the first incident shock wave, 
between the first reflected and second incident shock waves, and downstream of 
the second reflected shock wave were calculated from the pitot profiles assum- 
ing isoenergetic flow. A calibrated venturi meter was used to measure the 
bleed flow. 


ANALYSIS 

Figure 2 shows the flow model used in the analysis, R is the radial dis- 
tance from the tunnel centerline and x is the distance downstream from the tip 
of the shock wave generator. Conditions at station 1 are assumed to be known. 
The object of the analysis, given the shock generator shape and initial condi- 
tions, is to compute the locations of the reflected shock waves and the boun- 
dary layer properties at successive stations along the wall. Reflected compres- 
sion waves are treated in the analysis as a single shock wave. 

The boundary layer can be divided into three subregions associated with 
the three steps in which the computation is carried out. 

(1) Region I extends from Station 1 where the first incident shock wave reaches 
the boundary layer edge to Station 3 where the reflected shock wave emerges 
from the boundary layer. Surface 2 is the stream surface which passes through 
the intersection of the reflected shock wave with the boundary layer edge. Note 
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that this stream surface intersects the incident shock wave outside the boun- 
dary layer edge (see Fig. 3). This choice of surface 2 allows for mass entrain- 
ment into the boundary layer through the interaction region. The location of 
this surface and the pressure distributions along it are obtained from an in- 
viscid conical flow solution. Using this surface, the tunnel wall, and the 
planes normal to the wall at 1 and 3 to define a control volume, a control vol- 
ume analysis of the region may be used to determine the length of the inter- 
action, and the boundary layer thickness and shape at 3. The method used here 
is quite similar to the methods used by Seebaugh [1] and Mathews [2] except 
for the velocity profile representations for the boundary layer. The velocity 
profile representation used at Stations 1 and 3 and throughout the analysis is 
the improved wall-wake velocity profile developed by Sun and Childs [3] (see 
Appendix A). For turbulent isoenergetic compressible boundary layer flow, 
the modified profile may be expressed in the form. 



sin {arc sin [1+ i ^ (in n + |- (1-n^)^ 

e 

In (l + (l-n®)^^)) - I ^ (1+ cosn^)]} 


where 

n_ 

K 



6U 

In (— ^) - 5.1 + 

V 

w 


0.614, 
a K ^ 


and a = 1 is assumed. 


( 1 ) 

( 2 ) 


One other variation on the earlier method has been incorporated into the 
present analysis. The flow direction in the boundary layer downstream of the 
interaction has been taken as the average of the value at the wall (i.e., 
zero) and that at y = 6o as determined from the inviscid solution. In the 
analyses by Seebaugh [1] and Mathews [2] the flow direction downstream of an 
interaction was taken to be parallel to the wall. Figures 3a and 3b show the 
control volumes used in the present analysis. Although boundary layer bleed 
was not employed at the first interaction site in the experimental studies, 
the control volumes shown do allow for that possibility. For the control 
volumes shown, the continuity equation may be expressed in the form. 


w 


V^r^E 


27rpuRdR 


w 


V*3 


2TTpuRdR + mg 


(3) 


while the x-momentum equation may be written as. 
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V'^r^E 


w 


2Trp^RdR - J 2Trp2RdR - J P 32 irRdR- 2 TrR^LT^ 


R„-«r‘E 


w 




w 


R -6o 
w 3 


2irpu2RdR - 


w 


V'^T'^e 


2irpu2RdR + I 


Bx 


( 4 ) 


where 


’w ' \ * >"'* P; 


is the average static pressure over the boundary layer at 3. Comparable equa- 
tions may be written for the second interaction site. Allowance for boundary 
layer mass entrainment in the interaction region is made by assuming an en- 
trainment rate equal to that for the flow upstream of an interaction. 

(2) Region II extends from Station 3 to Station 5. This is the region of 
boundary layer flow between the first reflected shock and the second incident 
shock. Since no shock interactions are present in this region, the axial pres- 
sure gradient is relatively low. Starting with conditions at 3 as determined 
by the control volume analysis of the first interaction, changes in the boun- 
dary layer properties to Station 5 are computed using a turbulent boundary 
layer program suggested by Paynter and Schuelle [4]. The program uses a wall- 
wake profile to represent the velocity profile and the entrainment function 
concept proposed by Green [5] to solve the boundary layer equations (see 
Appendix B). An inviscid flow solution is used to provide the wall static pres- 
sure distribution needed for the boundary layer solution in this region. The 
inviscid solution, however, is obtained in a manner which allows for the 
effects of the first shock interaction. The method employed was to use an 
artificial wall position in the interaction region which would cause the re- 
flected inviscid shock wave position to match that determined by control volume 
analysis of region I. The wall static pressure distribution p(x) as determined 
using the artificial wall position is then assumed to exist at the actual wall 
and is used in the boundary layer calculations. 

(3) Region III extends from Station 5 to Station 7. This region covers the 
interaction of the second incident shock wave with the boundary layer. The 
method here is similar to that for Region I, except that the flow at control 
surface 6 is no longer conical. Conditions along this surface were obtained 
from a method of characteristics solution for flow past the double-cone center- 
body. In the characteristics solution the interaction of the first reflected 
shock with the second incident shock must be considered. The location in the 
flow field of the reflected shock wave was determined using the artificial wall 
position. 
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In the control volume equations (3) and (4) the mass bleed rate in mg is 
assumed to be known. The x-momentum of the bleed flow depends on the manner 
in which the bleed flow is accomplished. In the analyses by Seebaugh [1] and 
Mathews [2] computations were made for three bleed models: porous wall suc- 

tion, slot suction and scoop suction. Figure 3a shows the porous wall model. 
With this model, the x-momentum of the bleed flow, Ibx» was assumed to be zero 
Figure 3b shows the slot suction model. With slot suction, ig^ was assumed to 
have the same value as that possessed by the bleed mass as it entered the con- 
trol volume, i.e.. 


I 


Bx 



2TTpu2RdR 


(5) 


where Rg 


is determined from. 



27rpuRdR 


( 6 ) 


For the Mach 2.82 flow with bleed, the bleed holes were drilled normal to 
the wind tunnel wall. At first thought, then, the porous wall model might 
appear to provide a better representation of the bleed flow. However, the 
bleed hole diameter of 0.132 cm was on the order of one-half the boundary thick- 
ness at Station 3. Thus, there should be x-momentum associated with the bleed 
flow and the bleed flow behavior might then be expected to be somewhere between 
that for porous-wall suction and slot suction. As will be discussed in the 
section on results, this appears to have been the case. 


DESCRIPTION OF COMPUTATIONAL PROCEDURE 

The computation of the flow at M<„ = 2.82 with 2.8 percent bleed rate will 
be used here to illustrate the computational procedure. The computation is 
carried out in the following steps: 

(1) The boundary layer thickness and velocity profile shape at Station 1 are 
known from measurements. For the example under consideration the first inci- 
dent shock wave is at an angle of 22.75 degrees to the axis of the tunnel and 
reaches the boundary layer edge at a station approximately 5.2 cm downstream of 
the tip of the shock wave generator. Measurements were not taken exactly at 
Station 1. However, measurements were taken just upstream at x = 5.08 cm. We 
then assume that the boundary layer properties at Station 1 can be approximated 
by those obtained at x = 5.08 cm. The best representation of the boundary layer 
velocity distribution by the wall-wake profile, Eq. (1), is obtained by using 
the computer program, LEAST, listed in TABLE 1-A. 

In program LEAST, the boundary layer thickness 6 and the coefficient of 
wall friction Cf are regarded as two parameters and a double-loop-iteration 
scheme is used to solve for a least squares fit of the wall-wake profile to the 
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experimental velocity distribution. This program has been written for either 
isoenergetic or non-isoenergetic flow. The velocity distribution can be in- 
put in the form of pitot tube pressure or Mach number distribution. The in- 
put data listed in TABLE 1-B are pitot pressure profile data. Isoenergetic 
flow has been assumed. The description of the input is detailed in the com- 
ment cards listed in the program LEAST. For the case under consideration, 
the following boundary layer properties were obtained: 

6 = .406 cm, = 0.002, u^/u| = 0.0445 

and the freestream Reynolds number Re = 19.0 x 10^ per meter. 

(2) The control surface 2 is defined to coincide with the stream surface which 
extends upstream from the point of intersection of the reflected shock wave 
with the boundary layer edge. In the inviscid flow field behind a conical 
shock wave, the characteristic surfaces are conical surfaces originated from 
the tip of the shock generator. Therefore, the characteristics of the flow 
field are functions only of the half angles of the conical surfaces. Gootzait 
[6] developed a computer program, program C0NE, which is used to compute the 
flow field characteristics. The program C0NE and the input for it are listed 
in TABLES 2-A and 2-B. There is only one input card for the program. The in- 
put FORMAT is 7F10.6 and the input data are, in order, the specific heat ratio 
of the gas, the freestream Mach number, the half angle of the cone and the dis- 
tance, in inches, behind the tip of the cone where the output is desired. 

The conical stream surface can be described in non-dimensional form as. 


r 



+ 


f^/"o tan g 
^ tan <j>^ 


d(x/x^) 


the reference point (x„, rp) can be any point on the shock wave as shown in 
Figure 4a. a is the flow deflection angle at the stream surface as measured 
with respect to the tunnel axis and (jg is the incident shock wave angle. (Note: 
the final reference radius is obtained by an iterative process which allows for 
boundary layer mass entrainment. This is taken care of automatically in the 
program). Figure 4b shows the choice of reference radius rQ and the conversion 
of the stream surface to control surface 2. In the numerical computation a 
finite number of points is used to represent the control surface and the flow 
conditions such as Mach number, pressure and flow angle on each point are ob- 
tained from program C0NE as functions of characteristic angle (t>. 

(3) The computer program ANAL for the control volume, originally developed 
by Seebaugh [7], modified by Mathews [2], and then improved upon by the pre- 
sent authors is used to compute boundary layer properties at Station 3. The 
program ANAL is listed in TABLE 3-A, the inputs in TABLE 3-B. The input values 
at Station 1 for u./u*, 6-, and the Reynolds number based on boundary layer 
thickness, Ue«]/ve’ ®are obtained in step 1 and the conditions on the control 
surface 2 are obtained from step 2. The flow is assumed to be isoenergetic; 
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therefore the recovery factor is assumed to be 1. There is no boundary layer 
suction in this region. The shear stress at the wall between stations 1 and 
3 is assumed to be the average of the value at stations 1 and 3. The flow 
direction at station 3 is taken as the average of the values at the wall 
and at y = 6^ as determined from the inviscid solution. 

In this computation, we use the shock wave angle <f> , instead of the 
flow deflection angle across the first incident shock. We also choose to 
input conical flow conditions behind the first incident shock wave. The com- 
puted results include the interaction length and the boundary layer proper- 
ties at station 3. The predicted boundary layer properties at station 3 
include the boundary layer thickness, the coefficient of skin friction, the 
average static pressure across the boundary layer and the velocity distribu- 
tion. 

(4) The program ANAL, which is based on the control volume method, predicts 
only the properties of the boundary layer at the downstream and of the inter- 
action region. For computing the boundary layer downstream of the interaction 
and for computations of the next interaction, the flow field outside the boun- 
dary layer must be known. This is obtained by a method of characteristics 
computer program. The program used here is a modified version of a program 
originally prepared by Cavalleri [8]. The program designates the region be- 
fore the incident shock wave as region 2, the region between the incident re- 
flected shock waves as region 3 and the region downstream of the reflected 
shock wave as region 4. The flow in region 2 is uniform and parallel to the 
centerline and the flow in region 3 is conical flow. Figure 5 shows the desig- 
nated regions and input points at initial station x = 1.835 inches (4.65 cm) 

as used in the example under consideration. In computing the external flow 
field allowance must be made for the effect of the interaction. Otherwise 
the location of the reflected shock wave will be downstream of the actual re- 
flected wave and the external flow field needed for downstream boundary layer 
calculations will be in error. In order to avoid this problem, an artificial 
wall location which would cause the reflected shock location to match that 
predicted by the control volume analysis was used. Determination of the 
aoDropriate artificial wall location is a trial and error process and for the 
computations which have been carried out to date has required 

three runs of the characteristics program. In this computation it was found 
that a tunnel radius of 0.99 inches (2.51 cm) caused the reflected shock wave 
position to match that determined by the control volume analysis. The inputs 
to the program are listed in Table 4. The reader should refer to reference 8 
for the details of the inputs. 

(5) From station 3 to station 5 we use a turbulent boundary layer program 
BLGRN (see Appendix B) which is a modified version of one developed earlier 
by Paynter and Schuelle [4] based on Green's [5] entrainment function concept. 
The control volume analysis in step 3 provides the initial conditions at 
station 3. For the boundary conditions on surface 4, we assume that the Mach 
number and pressure distributions are equal to those obtained along the arti- 
ficial wall in the inviscid solution of step 4. We should also point out here 
that the computation is a two-dimensional approximation. However, in view of 
the ratio of boundary layer thickness to tunnel radius the results should be 


9 


applicable for the ax i symmetric flow under consideration. The program, BLGRN, 
is listed in TABLE 5-A, the inputs in TABLE 5-B. The inputs are described 
in the comment cards at the beginning of the program listing. 

(6) At 1.942 inches (4.94 cm) behind the tip, the shock generator broke from 
a 10 degree cone to 13 degrees and generated a second shock wave. The posi- 
tion of the second incident shock wave and the flow field characteristics be- 
hind the second shock wave are computed by using the same method of character- 
istics program which is used in step 4. Given the geometry of the centerbody 
and the inviscid solution (based on the artificial wall position) from step 4, 
we have the inputs for the computer program as listed in TABLE 6 and illustrat- 
ed in Figure 6. As in the computation of step 4, the program designates the 
region before the second incident shock wave as region 2, the region between 
the second incident and reflected shock waves as region 3, and the region 
downstream of the second reflected shock wave as region 4. In region 2, un- 
like the situation in step 4 where the external flow is uniform, there exists 

a reflected shock wave. The program is not capable of solving for the inter- 
section of two shock waves. In preparing for the input, the first reflected 
shock wave is replaced by a smoothed compression wave band. In the calcula- 
tion described here a linear compression spaced over six input points was 
used. The flow between the artificial wall radius of 0.99 inch (2.51 cm) and 
the actual wall radius of 1.02 inches (2.60 cm) was assumed to be uniform and para- 
lel to the wall. The result shows that the second incident shock wave has an 
angle of approximately 33 degrees and that it reaches the boundary layer edge 
at x = 2.85 inches (7.22 cm) where the boundary layer thickness is 0.125 
(0.312 cm) as computed in step 5. 

(7) The control surface 6 is also defined to lie along a stream surface be- 
hind the second incident shock wave, but the flow field in this region is quite 
complicated. Because of the nature of the characteristics program, it was 
necessary to run the program three separate times to determine flow properties 
at the three stations. The pressure and Mach number distributions at each sta- 
tion are obtained from the result of the previous step. At x = 2.85 inches 
(7.22 cm) where the second incident shock wave reaches the edge of boundary 
layer the cone radius is 0.552 inch (1.40 cm) while the distance from the cen- 
terline to the edge of the boundary layer is 0.895 inch (2.28 cm). The mass 
flux through the area between the radii is 0.412 Ibm/sec (0.187 Kg/sec.). 

At X = 2.988 and 3.205 inches (7.60 and 8.15 cm) for the mass flux to be equal 
to the 0.412 Ibm/sec (0.187 Kg/sec) the radii are 0.907 and 0.922 inch (2.30 
and 2.34 cm) respectively. Ideally, we can compute at more stations to de- 
termine the stream surface. But the preparation of pressure and Mach number 
profiles is very tedious. We determine the stream surface by graphically 
tracing a smooth curve through the above three points. The computer program 
MFLX is used here to compute mass flux from the profile of pressure and Mach 
numbers. Program MFLX and the input to the program are listed in TABLE 7. 

(8) The computer program ANAL listed in TABLE 3-A is used for computation 
from station 5 to station 7. The boundary layer properties at station 5 are 
obtained from step 5 and the location of the control surface 6 is determined 
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by, step 7. After determining the control surface 6, the flow properties along 
the surface are obtained from the output of step 6. The x-coordinate should 
be transformed into the distance to the axial station at which the extension 
of the second shock wave intersects with the tunnel centerline as shown in 
Figure 8. For the points on the control surface the input to the program, 
i.e., the x-coordinates, should be normalized by distance Xq and the r-coordi- 
nates should be normalized by radius Rq. The distance Xq is the distance from 
station 5 to the intersection point of the second shock wave with the center- 
line and the radius Rq is the distance from the tunnel centerline to the edge 
of the boundary layer at station 5. For the flow considered in the sample cal- 
culations boundary layer bleed takes place at the second interaction site. The 
bleed mass flux is about 2.8 percent of the boundary layer mass flux at station 
5. In this computation the slot suction model has been assumed for the bleed 
flow. The input to the program ANAL is listed in TABLE 8. The details of the 
input are described in the comment cards at the beginning of program ANAL list- 
ed in TABLE 3-A. 


RESULTS 

The results of the analysis for three different flow cases, one at = 3.82 
with no bleed, two at = 2.82 with bleed at the second interaction site, are shown 
in Figures 9 through 16? Results from the sample calculation described herein are 
given in Figures 12 and 13. Comparisons are made between predicted and measured 
results. The data for M = 2.82 flow are from an investigation by Teeter [9]. 

Figure 9 shows comparisons of the experimental and predicted shock wave 
patterns and boundary layer thickness for the Mach 2.82 flow. Since the pre- 
dicted shock wave locations are determined by the inviscid analysis used in 
combination with the artificial wall position, no induced shock wave is pre- 
dicted. Also shown is the pressure distribution at the tunnel side wall as a 
function of the distance aft of the cone tip. The triangular points shown for 
the analysis in the static pressure plot were determined by using the artifi- 
cial wall and the inviscid flow solution. The predicted and observed values 
are seen to be in good agreement along the entire length of the double shock 
interaction. 

* 

Figure 10 shows 6 , e and at several stations along the tunnel side 
v/all. Here also, predicted and observed results are in good agreement. The 
experimental values shown for Cf have been determined by a least-squares fit 
of the modified wall-wake velocity profile to the experimentally determined 
velocity profiles. 

Figure 11 shows Mach number profiles for the boundary layer at the down- 
stream end of the second interaction. The analysis predicts the end of the 
second interaction to be at x = 3.75 inches. Since profiles were not taken at 
this specific station, profiles taken just upstream at x = 3.70 inches (9.5 cm) 
and just downstream at x = 3.80 inches (9.65 cm) are shown for comparison. It 
is apparent that the analysis leads to a profile which provides a good repre- 
sentation of the experimentally determined profile near the interaction end. 
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Figure 12 shows boundary-layer thicknesses, shock wave patterns and wall 
static pressure distributions for the Mach 2.82 flow with 2.8 percent boun- 
dary-layer bleed at the second interaction site. Figure 13 shows 6*, e, and 
for this flow. 

Two sets of predicted results are given, one for the porous-wall suc- 
tion model, the other for slot suction. (For the sample calculation described 
in the report slot suction has been assumed). As is shown, the differences 
between the results for the two suction models are not large. Differences in 
predicted results with the two models are due solely to the differences in 
values assigned to the x-momentum of the bleed flux. Since the bleed rate is 
lew, the x-momentum associated with the slot suction model is small and not too 
different from the zero values for porous suction. The predicted and measured 
results are in reasonably good agreement. 

Figure 14 shows boundary-layer thickness, shock wave patterns and wall 
static pressure distributions for M = 2.82 with 5.0 percent bleed. Values for 
6*, e, and Cf are shown in Figure 15, while Mach number profiles downstream 
of the second interaction are shown in Figure 16. The Mach number profiles 
represented by the solid lines in Figure 16 are predicted profiles for the two 
bleed flow models. They are shown on the figure at the axial positions pre- 
dicted for the end of the second interaction. Since experimental profiles 
were not taken at these precise locations, experimental profiles taken in the 
neighborhood (x = 3.20, 3.30 and 3.40 inches or 8.13, 8.38 and 8.64 cm.) of the 
predicted locations have been shown for comparison. 

The flow conditions up to the second interaction are the same as those 
for the flow with 2.8 percent bleed. With the higher bleed rate the differ- 
ence between the results for porous wall and slot suction are much more pro- 
nounced than with 2.8 percent. The slot-suction model gives a reflected shock 
location which is in better agreement with the observed results. On the other 
hand, the values of 6*, e, and Cf obtained with the porous wall model agree 
better with experimental values than do the slot suction results. As was 
pointed out in the section on analysis, the bleed hole diameter of 0.132 cm 
was on the order of one-half of the boundary-layer thickness so that the 
bleed flow behavior might be expected to lie between that for porous wall 
suction and slot suction. The x-momentum of the bleed value might then, in 
turn, be expected to lie between the values used with the two models. Indeed, 
the use of a bleed flow momentum flux between the two limits would lead to 
better overall agreement between predicted and measured values of 6*, e, and 
Cf. Even then, however, the predicted interaction length would be too long. 

It should be remarked that in estimating Ibx foi" the slot-suction bleed flow 
model, no allowance is made for the turbulent shear stress along the stream 
surface separating the bleed flow from the main body of the flow, nor for the 
wall shear stress. Nor is the pressure force along the separating stream 
surface considered. The effects of the pressure force and wall shear tend 
to cancel the effect of the turbulent shear on the separating stream surface, 
but the extent to which they do so is not known. It should be remarked fur- 
ther that no allowance is made for the roughness effect of the holes on the 
wall shear. Further study is needed on the details of the bleed flow behavior, 
including the roughness effect of the holes, before the effects of bleed con- 
figuration can be resolved. 
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The computation procedure reported here represents the results of a con- 
tinuing effort to improve analytical methods of predicting flowfields in the 
inlet of supersonic aircraft. In a recent analysis of inlet flowfields by 
Reyhner and Hickcox [10] the effect of the shock wave interaction on the in- 
viscid flow was taken into account by first obtaining a control volume solu- 
tion for the boundary-layer properties downstream of the interaction. Then, 
using an effective surface defined by the boundary-layer displacement thick- 
ness upstream and downstream of the interaction, and using a patching techni- 
que across the interaction region to construct an effective displacement sur- 
face for that region, the inviscid flow solution was obtained for the effective 
surface. A comparable technique was tried in the work reported here but it 
was not as successful as the scheme of using the simple reflection off the 
artificial wall. 


CONCLUSIONS 

A control volume analysis method, employed in conjunction with a turbu- 
lent boundary-layer computation scheme, has been used to predict the flow- 
field downstream of successive shock wave boundary-layer interactions for 
flows at Mou = 3.82 and 2.82. The computational procedure has been outlined 
in detail. The effects of boundary layer bleed at the second interaction site 
have been considered. For flow with low bleed rates or no bleed the predicted 
interaction lengths and wall static pressures, as well as the boundary-layer 
properties downstream of the interactions show good agreement with measured 
results. With low bleed flow the predicted results for the slot-suction and 
porous-wall models differ only slightly since the momentum of the bleed flow 
is small. As the bleed flow rate is increased, predicted and measured re- 
sults are also in reasonably good agreement. Here, however, differences be- 
tween predicted results for the two suction models are larger since the dif- 
ference between the momentum fluxes of the bleed flows is larger. A value of 
bleed flow momentum between the values used for the models would improve the 
agreement between predicted and measured results. 
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APPENDIX A 


Modified Wall -Wake Profile 

A simple representation of the mean velocity distribution in a turbulent 
boundary layer is very useful in integral analyses of turbulent flow problems. 
After an extensive survey of mean velocity profile measurement, Coles [11] sug- 
gested that for incompressible turbulent boundary layer flow the velocity pro- 
file may be represented by a linear combination of two universal functions in 


the form. 

u/u 

T 

= (1/K)ln (yu^/6)+C+n W (y/6)/K = f(y)+g(y) 

(A-1) 

where 

f(y) 

= (1/K)ln (yu^6) + C 

(A-2) 

is the Law 

of the Wall and 


g(y) 

= nW (y/6)/K 

(A-3) 


is the Law of the Wake. 

Setting u/Ug and W(y/6) = 2 at (y/6) = 1 in Eq. (A-1) and subtracting 
the resulting equation from Eq. (A-1) leads to an expression for the velocity 
of the form, 

u/u^ = 1+(1/K)(u /uJln(y/6) - (n/K)(u /uJ[(2-W(y/6)] (A-4) 

c T c T c 

Mathews et al., [12] have developed a wall-wake representation of the 
velocity profile in a form applicable for isoenergetic compressible boundary 
layers. Their profile is expressed as, 

u/Ug = (l/o'^)sin {arc sin [l + (l/K)(u^/u*)ln (y/6) 

- (n/K)(u^/u*)(l + cos (iry/6))]} (A-5) 

where 

(u^/u*) = [(Cf/2) a/(l-a)]^Varc sin (A-6) 

and 

n/K = (1) {[l/(u^/u*)] - (1/K)ln [Re^ (Cf/2)^^(l-a)^ *2®] - C} (A-7) 

and where (2-W) has been replaced by 1 + cos (iry/6) for mathematical con- 
venience. 
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Equation (A-5) has been found to provide a good representation of the 
boundary layer velocity profile for a range of external Mach numbers and wall 
static pressure gradients. However, with both Eq. (A-4) and Eq. (A-5) the 
velocity gradient at the boundary layer edge is found to have a non-zero 
value. In the modified wall -wake which is to be developed here this short- 
coming is avoided. 

The law of the wall may be derived from Prandtl's mixing length theory 
and the assumption that the shear stress is constant across the boundary layer 


(Cf. Schlichting [13]). However, an expression for t of the form, 

T = [1 - (y/6)®] = (A-8) 

Where a is a real constant should provide a more realistic relationship 
for the shear stress. 

Using Eq. (A-8) we may write, 

Tw (1-n®) = pK^y2 (du/dy)2 (A-9) 

Integration of Eq. (A-9) gives an expression for u/u^ of the form, 

u/u^ = (1/K)ln n + (2/aK) {(1-n®)^ - In [l+(l-n®)^]} + (A-10) 

Replacing f(y) in Eq. (A-1) by Eq. (A-10) we have, 

u/u^ = (1/K)ln n + (2/aK) {(1-n®)*'" - In [1+ (1-n®)^]} + 

+ (n/K) W(n) (A-11) 

At the boundary layer edge (n^l) we have, 

Ue/u^ = Cl + (n/K) W(l) = + 2 n/K (A-12) 

while near the wall (as n-^0^). 


u/u^ = (1/K) In (yu^/v) - (1/K) In (6u^/v) + + 0.614/aK (A-13) 

Near the wall the expression for the law of the wall as given by Eq. 
(A-2) is also applicable. Equating the expression for u/u we may evaluate 


= 5.1 - (0.614/aK) + (1/K) In (6u^/v) 


(A-14) 
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while from Eq. (A-12) 

n/K = (1/2) [(Ug/u^) - 5.1 - (1/K) In (6u^/v) + 0.614/aK] (A-15) 

Following procedures similar to those used by Van Driest [14], Maise 
and McDonald [15] or Mathews [12], 


u. 


(B2+4A2)' 

~2F 


sin {arc sin 


2A2-B n . 1 % ,, 

TT n + F 77* ' ^ 


(B2+4A2)' 


K u* 
e 


+ ' - I In (l + (l-n^)^^) - f ^ (2-W(n))]} + 

e 


B 

W 


(A-16) 


where 


n/K 


(1/2) [(u*/u^) - (1/K) In (6U^/V^) - 5.1 + 0.614/aK] 


(A-17) 


and 


u*/u^ = [(2A2-B)/(B2+4A2)*^] 


(A-18) 


For isoenergetic flow, Eqs. (A-16) and (A-18) become, respectively. 


and 


— = -V sin (arc sin 

e o 


n+ ^ ^ (Inn + 


2(l-n^) 


9 


- flnd + d-n^)"'^)) 


K 



( 1 + cos Tin) ] } 


(A-19) 


u /u* = [(Cf/2) a/(l - o)]'^/arc sin a'^ 

T G T 


(A-20) 


where 2-W(n) has been replaced by 1 + cos nn for mathematical convenience. 
As a •* «> Eq. (A-11) reduces to Eq. (A-1) while Eq. (A-19) reduces to the 
profile proposed by Mathews et al . , Eq. (A-5). 


The remaining problem is the selection of the constant a. Based on mea- 
surements reported by Klebanoff [16] and Horstman and Owen [17], it appears 
that a = 1 represents a reasonable assumption. This amounts to the assumption 
of a linear shear stress distribution across the boundary layer. 


The method of least squares has been used to fit the wall -wake profile, 
Eq. (A-19), for both a = 1 and a ^ “ to a number of experimental velocity 
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profiles by Seebaugh [7], Teeter [9] and Rose [18]. The computations can be 
carried out by using program LEAST listed in TABLE 1-A. An example of the 
results is given in Figure 17 which shows two profiles from the study by 
Seebaugh of an interaction between a conical shock wave and the turbulent 
boundary layer at the wall of an axially symmetric M = 2.82 wind tunnel. The 
profiles are for stations just upstream and just downstream of the interaction 
region and the experimental velocities have been calculated from pitot pres- 
sure and the wall static pressures under the assumption of isoenergetic flow. 
The 10-degree half angle cone used in the study did not produce a shock wave 
of sufficient strength to cause boundary layer separation. The values of Cf 
and 6 determined by the curve fits are listed on the figure along with values 
for the displacement and momentum thickness 6* and e. 

As is shown in the figure, both the modified wall-wake profile (a=l) and 
the profile for a ^ ~ provide good representations of the experimental velo- 
city distribution over the ranges from y = 0 to the values determined for 6 by 
the curve fits. The values of Cf, 6* and e determined for the two profiles 
differ only slightly. However, the values of 6 as determined with the modified 
profile show much better agreement with the values of 6 based on u/Ug = 0.995. 
Furthermore, the velocity gradient for the modified profile goes to zero at 
y = 6. In all of the data examined to date the modified velocity distribution 
than the earlier version of the compressible wall-wake profile. 


APPENDIX B 


Boundary Layer Computations 

For turbulent boundary layers in compressible flow. Green [5] derived 
from Head's (Cf. Green [5] ) work for incompressible flow a procedure for simu 
lataneously calculating the development of the momentum thickness and a quan- 
tity referred to as mass flow thickness, defined as, 


A 


6 

‘ 0 


PU 


P u 
^e e 


dy 


6 - 6 * 


(B-1) 


In this procedure the momentum-integral equation, 
Cx 

e' 


dx 


5 a du^ 

■ r e d3T 


. e dr 


r dx 


(B-2) 


is integrated simultaneously with an auxiliary equation which accounts for the 
rate at which the boundary layer entrains fluid from the free stream. 


d3T 


e e 




(B-3) 


In equation (B-2) the value of j is set equal to zero for two-dimension- 
al flow, to unity for axi symmetric flow. 


Equation (B-3) can be rearranged in the form. 


f. - ^ * K 


A du 

1 ) — — ^ 
' u^ dx 
e 


(B-4) 


Green [5] found semi -empiri cal ly that F has the following form, 

,-0.653 


0.0306 ((H-|),< - 3.0)' 


(B-5) 


where 


(H,) 


^dy 

'0 ^e 


VK 


6 U 


^ 0 - ^) dy 
e 


(B-6) 
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Paynter and Schuehle [4], using empirical expressions for (Hi)|/ and Cf 
suggested by Green, developed a computer program to solve equations (B-2), 
(B-4) and (B-5), However, as it is suggested by Sun and Childs [3] the velo- 
city in compressible turbulent boundary layer flow can be represented by equa- 
tion (1). Using equation (1), we can solve equations (B-2), (B-4) and (B-5) 
without dependence on empirical formulas for (H-|)|^ and C^, 


From equation (B-1) we have, 

dC. 


dx 


dA 

dx 


36 

3C, 


dx 


39 

36 


dx 


dM 


36* ^ 

3Ci dx 


f36* 

^36 


39 '"‘e 

3M„ dx 
e 


(B-7) 

d6 36* 

dx 3M^ 

e 

dM„ 

e 

dx 

(B-8) 


dC- 


Solving equations (B-7) and (B-8) for and ^ yields. 


dx 


dC^ 

dF 


, , 36* S / de 36 

' 36 '^dx ■ 3M 


dM 
€ 

dx 


/ 39 \ / dA 
^16^ ^dx 


36* 

3M^ dx ‘ 
e 


30 ^ 96* ^ 

^ ‘ ‘ 36 ^ 3C^ 9C^ 36 


(B-9) 


dx 


90 


rdA 


3C^ ^dx 


36* *^^e ^ . 3^ _ 3^ ^ 

3M^ dx ' 30^: ^dx " 3M. dx 

e t e 


n I^) 99 + 36* 39 
' ■ 36 ^ 3C^ 3C^ 36 


(B-10) 


To evaluate the displacement thickness and momentum thickness and their 
partial derivatives with respect to 6, Cf and Mg we use the expression given 
by Van Driest [14] for the temperature distribution through the boundary 
layer. 


— = — + (1 + _ Zw) i!_ _ 111 (— )2 

e e e e e 

Equation (B-11) was used in the derivation of equation (1). 


(B-11) 


The displacement thickness and its partial derivatives with respect to 
6, C^ and Mg can be written as. 


6 * 


(1 . ^) dy = 


"e"e 


- TO-) 


(B-12) 
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36* 

3C. 


36* 

3M„ 


(1 - T ^^) I 

e y=5 


— (1 

36 




^ 9_ fZ§. iL_^ riv 
0 e 


^^3 n u X H./ 

^ (' " T"T7"^ 

0 “e 


(t^ ^■) dy 


0 ^ 


3 / e u \ 

0 ’«e ‘t u,) 


(B-13) 


(B-14) 


(B-15) 


Similarly, we have the momentum thickness and its partial derivatives 
as. 
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PU 


0 Pe^'e 


(1 - dy 


T u 

^ 0 - ^) “y 


30 

36 




■) dy - 


•6 

0 


3 / ^e 

36 Ug 


) 


dy 



j-6 

^0 




dy - 


'6 

-0 



e iJ 
T u 


r) 


dy 





dy 



dy 


(B-16) 


(B-17) 


(B-18) 


(B-19) 


If initial values of <5, Cf and Mg are known an initial velocity profile 
can be determined by substituting values of '5, Cf and Mg into equation (1) and 
values of ® and and their partial derivatives can be determined. Values of 
(Hi)k can be computed by equation (B-6). We then use the values obtained and 
equations (B-2), (B-4) and (B-5) to compute d /dx and d /dx. Substituting the 
obtained values into equations (B-9) and (B-10) gives values for dCf/dx and 
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I 


dCf/dx. The computations for s and Cf may be carried out in step-by-step fashion. 
Let Ax be a small increment of x; then values of s and Cf at x + ax may be 
expressed as. 


6 (x + Ax) = 

6(X) + ^ AX 

(B-20) 

C^ (X + Ax) = 

dC, 

(B-21) 


The boundary conditions Mg and dMg/dx are assumed to be known, or may be 
determined from the method of characteristics solution. By repeating the com- 
putation, the properties can be computed throughout the boundary layer flow. 
Program BL6RN, listed in TABLE 5-A, has been developed for the numerical com- 
putations. A sample of the results computed with this program is shown in 
Hirst [19]. As is shown, derivations between computed and experimental values 
of H, and Ree begin to show up as a condition of separation is approached. 
However, the computed values of Cf agree quite well with the data over the 
entire range of the computations. 
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(a) POROUS WALL BLEED 



(b) SLOT BLEED 


FIGURE 3. CONTROL VOLUME USED IN ANALYSIS 
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X INPUT DATA POINTS TO METHOD 
OF CHARACTERISTIC PROGRAM 



FIGURE 5. INPUT FOR COMPUTING EXTERNAL FLOW FIELD 
BY METHOD OF CHARACTERISTICS 


y INPUT POINTS. TO METHOD OF 

characteristic program 



FIGURE 6. INPUT FOR COMPUTING EXTERNAL FLOW FIELD 
FOR SECOND INTERACTION 
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X 1NPBT DATA JOINTS 
TO FROGNAH MFLX 



FIGURE 7. COMPUTATION OF MASS FLUX AS A FUNCTION 
OF RADIAL COORDINATE 


X INPUT DATA POINTS FOR BOUNDARY 
LAYER EDGE STREAM LINE 



FIGURE 8. LOCATION OF STREAM SURFACES 
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FIGURE 9. SHOCK WAVE POSITIONS, BOUNDARY LAYER THICKNESS 
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FIGURE 14. SHOCK WAVE POSITIONS. BOUNDARY LAYER THICKNESS 
AND WALL STATIC PRESSURES. M = 2.82, 5.0 
PERCENT BLEED. 
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FIGURE 17. VELOCITY PROFILES UPSTREAM AND DOWNSTREAM OF 
A SHOCK WAVE - BOUNDARY LAYER INTERACTION. 
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TABLE 1-A 


PROGRAM LEAST 


PROGRAM LEAST (INPUT, OUTPUT, TAPE 5=INPUT, TAPE 6=0UTPUT) 


PROGRAM FOR LEAST SQUARES FIT OF WALL HAKE PROFILE 
INPUT FORMAT 7F13.6 EXCEPT CARD 1 
CAROCS) COLUMNS 

1 TITLE, COLUMNS 1-72 HOLLERITH 

2 1-10 AT =0. NO TOTAL TEMPERATURE DISTRIBUTION INPUT 

=1. total TEMPERATURE DISTRIBUTION INPUT 
11-20 AIP =0. CONSTANT PRESSURE DISTRIBUTION 
=1. LINEAR PRESSURE DISTRIBUTION 
=2. POINTHISE PRESSURE INPUT 
21-30 APT =0. NO PITOT PRESSURE INPUT 

=1. PITOT PRESSURE DISTRIBUTION INPUT 
31-40 AIM =0. NO MACH NUMBER DISTRIBUTION INPUT 
=1. MACH NUMBER DISTRIBUTION INPUT 
41-50 AIY =0. EQUAL Y-INCREMENT 

=1. POINTHISE Y INPUT 

51-60 AJB =0. NO MORE JOB AFTER THIS INPUT 
=1. MORE JOB AFTER THIS INPUT 

3 1-10 CA =5.1 

11-20 BK =0.4 

21-30 AID =1. 

31-40 AAK =5C. 

41-50 UTFT =0.04 FIRST ASSUMED VALUE FOR UT 

4 1-10 GAMMA =1.4 GAS CONSTANT 

11-20 PTOl TOTAL PRESSURE OUTSIDE B.L. IN PSIA 

21-30 TTOl TOTAL TEMPERATURE OUTSIDE B.L. IN R 

5 1-10 AM NUMBER OF PROFILES INPUT 

6 1-10 PH STATIC PRESSURE AT WALL 

11-20 R LOCAL RADIUS OF TUNNEL IN INCH 

21-30 X STATION DISTANCE 

31-40 AN NUMBER OF POINT INPUT 

41-50 ANN APPROXIMATE NUMBER OF POINTS TO BL EDGE 
51-60 OY SIZE OF Y-INCREMENT IN INCH 

7 1-70 PT PITOT PRESSURE.AN VALUES IN PSIA (APT=1.) 

8 1-70 TT TOTAL TEMPERATURE AN VALUES IN R (AT=1.} 

9 1-70 EM MACH NUMBER AN VALUES (AIM=1.) 

10 1-70 Y Y INPUT IN I INCH, AN VALUE (AIY=1.I 

11 1-70 P STATIC PRESSURE INPUT (AIP=2.» 

DIMENSION PW(9> ,PE£ (9) ,PTU( 60,9 J,PTU£ (9) ,PT (60,9) ,P (60, , 

1 TT(60,9),T(6C,9),TT£(9),EM(6C,9) ,RHO(60,9) , 

2 R(9> ,X( 9) ,ANN (9) ,Y(6C ,9) ,22(60) 

DIMENSION 2(60 ,UZ(60) , YZ(60) 

DIMENSION TITLE (12) , Y¥(60,9) 
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DIHENSION ZP(6D)tZMC601»ZTC60> 

DIMENSION ZS(60)»ZH(60> 

COMMON A1C*3K,CA«GAMMA, GAMl « G AH 2 , GAH3 , G AH<t t G AM5 « G AM6rEEH« 

1 ANC9)»CFF<9I,DEL(9),RH0EC9)»UE(9I,TEC9I» 

2 U(6t?,9) ,AAK,UTFr 

90 FEAOCEtlOlO) I TITLE ( J ) » J=l* 121 

R€ADC5»1300) AT, AlP » APT ♦ AIM » A lY, A JB 
READC5,lfl00) CA,BK,AIO,AAK,UTFT 
R£AO«5,10C0l GAMMA, PTOltTTOl 
GAM1=GAMMA-1, 

GAM2=GAMMA+1. 

GAH3=GAMMA/GAM1 

GAM4=1./GAM3 

GAM5=1./GAMMA 

6AM6=GAM1/2, 

GAM7=1,/GAM1 
READ(5,i000) AM 
M=AM 

DO 8 J=1,M 

1 R£AD(5,iOOOl PW(J1,R(JI,X(J1,ANCJ),ANN(J1,DY 
AN1=AN(J) 

AN2=ANN( J1 
N=AN1 
N2=AN2 
L=1 

IFCAPT.EQ.l.l REAO(5,1COO) (PTCI,J),I:=1,N 1 
IF(AT*EQ.l.) READ<5,1C0C) « TT ( I , Jl , 1=1, N ) 

IFfAIM.EO.l.l READI5,1C00J (EM (I , J ) , 1=1 , N1 
IF(AT,EQ,1.) GO TO 103 
DO 102 1=1, N 
TT(I,J1=TT01 
1C2 CONTINUE 

103 IFCAIY.EQ.O.) GO TO 104 
READ (5, 100 0) (Yd, Jl ,I=1,N) 

GO TO 223 

104 Y(1,J)=0. 

DO 2 1=2, N 

Y (I,J)=Y(I-l,J)*OY 

2 CONTINUE 

223 IF(AIP.EQ.2.) REAO(5,100Q) (P f I , J ) ,1=1 , N1 
0EL(J1=Y(N2,J1 
11 DO 3 1=1, N 

YY(I,JI=Y(I,J1/0ELCJ) 

Z(I) = YCI, Jl 
ZZ(I)=PT(I,J| 

ZP(II=P(I,J) 

ZT(I)=TTCI, J) 

ZM(I1=EM(I, J) 

3 CONTINUE 
DLL=DEL( Jl 

IFIAIM.EO.l.l GO TO 14 
IF(AIP.EQ.2.) GO TO 15 

COMPUTE EME BY PT02/PT01 

IF(L.EQ,1) PT02=PT(N2,JI 

IF(L.GT.l) CA4.L INIP(Z,ZZ,DLL,PT02,N) 
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PT2PT1=PT02/PT01 

EMOLO=1.666*PT2PTl»PT2PTi-5^666*PT2PTl*5. 

4 ARl=(GAM2*EM0La*EM0LD/C G AN l*EN0LD»EM0LDi-2 . ) ) **GAHMA 
AR2=PT2Pri»»GAMl 

AR3= ( AR1*GAM2/AR2+GAM1) /C 2.*GAMMA> 

EMN£W=SQRT(AR3I 

IF(ABS(ENOLO-EMNEWI-0.0 01> It). 10* 20 
20 EM0L0=EMNEM 
GO TO 4 

14 IFIL^EQ.l) EMNEM=EHtN2. J) 

IFtL.GT.l) CALL I NTP CZ. ZM.OLL.EHNEW.NI 

10 EH£=EMNEW 

COrtPUTE PE BY PTOl AND £N£ 

AR4= ( l.+GAMl*EME*EME/2. ) ♦♦GAM3 

PE=PT01/AR4 

PEE<JI=PE 

IF! AIP.EQ.2.1 GO TO 15 
DO 5 1=1. N 

P(I.J)=(P£E<JI-PWCJ) )»YYtI.J)*AIP^PWIJI 

5 CONTINUE 

IFCAIP.EQ.O.) P£=PWIJ) 

15 IFfAIP.EQ.2.) CALL INTP (Z .ZP. OLL.PE. Nl 
PEE<J)=PE 

TS=TT01 

IF(AT.Ea^l-) CALL INTP <Z .ZT. DLL. TS.NI 
TTE(J)=TS 

COMPUTE EM BY PT2/P 

IFIAIM.EQ^l.l GO TO 17 
TEST=(0.5*GAM2)»»GAM3 
00 6 1=1. N 

PATIO=PT(I.J)/P(I.JI 
IFIRATIO.LT.I.) GO TO 30 
IF(RATIO-TEST) 4C.50.60 

SUBSONIC 

40 EM(I. J)=SQRT(2.^IRATI0**GAN4-1* l/GAMl.) 

GO TO 6 

PT2 LESS THAN P 

30 EMfI.J)=0« 

GO TO 6 

SONIC 

50 EM<I.J)=1« 

GO TO 6 

SUPERSONIC 

60 EHOLO=1,0 5 
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63 A«l= < <2.*GAMM A*EM0L0»EM0LD-GAM1 )/GAH2) ♦♦GAMS 
AR2=RATI0^»GAM4 
tHNEW=SQRT(2.»ARl^AR2/GAM2l 
IFCABS(EMOLO-EMN£WI-0.001I 61f61»62 

61 £M(I,J)=EHNEM 
GO TO 6 

62 EMOLO=EMNEM 
GO TO 63 

6 CONTINUE 

COMPUTE UPSTREAM PROPERTIES 
IT DO 7 1=1, N 

PTU<I,JI=P(I,J)»(1.+GAM6^£M(I,J)»£M(I, JII*»GAM3 
TCI, J)=TT (I,J1/C1*+GAM6^EM(I, J) ♦EM(I,J1 ) 

RT=SQRTCT II, J) ) 

UII, J)=EM (I, J)*49.0»RT 

RHO(I,J) = F(I, J)^144 •♦32.2/C 1716. ♦TCI, Jll 

7 CONTINUE 

IFCAPT.EQ.l.l GO TO 177 
DO 117 1=1, N 
PTII, J)=PTUII,J) 

lEIEMlI, J) .GT.l.) PTII, JI=PII,J)^ CC.5»GAM2»EMfI,J)^EMII, J) I 
1»*GAM3/I2.*GAMMA»EM (I ,J)^EM II ,J ) /GAM 2-GAM 1/GAM2 ) ♦♦GAM7 
117 CONTINUE 

177 IFlAIM.EQ.l.) GO TO 16 

COMPUTE EDGE CONDITION 

TEST=IC.5^GAM2)^^GAM3 
PATI0=PT02/P££I Jl 
IFIRATIO.lt. l.i GO TO 39 
IFIRATIO-TESTI 49,59,69 

SUBSONIC 

49 EME =SQRTC2.»( RATI0^^GAM4-1. )/GAMl) 

GO TO 96 

PT2 LESS THAN P 

39 EME=0, 

GO TO 96 

SONIC 

59 EME=1. 

GO TO 96 

SUPERSONIC 

69 EMOLD=1.05 

163 AR1= I I2.»GAMMA»EM0LD»£M0L0-GAM1 )/GAN2l ♦♦GAM5 
AR2=RATI0^*GAM4 
EMNEW=S0RT(2.*ARi»AR2/GAM2J 
IFIABSIEMCLD-EMNEWl-O.Oui) 161,161,162 
162 EMOLD=EMNEW 


o o 


GO TO 163 
161 EHEsEHNEM 
96 CONTINUE 

COMPUTE UPSTREAM PROPERTIES 
16 PTUE( J)=PEE(J) »(1.+GAM6*EME*EMEI**GAM3 
TE(J)=TTE (J)/(1.+GAM6*EHE*EME> 

PT=SQRT(T£<JI > 

U£<J) =EME*49.*RT 

RHO£( JI=PEE(J) *144.*32.2/<1716.*TE<jn 
EEM=EME 

THTO=TT(l,J)/TT£(J) 

CALL LHLH (JtZtSNNtNOtLtTMTO) 

IP<ABS(SNN) .LE.O. OOOOOJOl) GO TO 21 

IF(L.GT.l) GO TO 12 

L=L+1 

DE1=0£LC J> 

0£L(J)=a. 9«OEL(J} 

SNN1=SNN 
GO TO 11 
12 L=L*1 

SLAPE=SNN1/<SNN-SNN1» 

0E2=DEL( J) 

OEL( J)=DEl*SLAPE*{OEl-OE2» 

D£1=DE2 
SNN1=SNN 
GO TO 11 

21 HRITE16,2010) (TITLE (I) t 1=1. 121 
HRITE(6,2003) X ( J) .PTOl.TTOl 
WRITE (6.2C04) 

DO 210 1=1. N 

WRITE (6.2 005) I.Y (I . J) . P ( I . J) .PT( I. J) . £M( I, J) . TT( I, J) 

210 CONTINUE 

WRITE(6.2010) (TITLE(I) .1=1.12) 

WRITE (6. 2 0 00) X(J) .R(J) .PW( J) .0£L(J) .CEP ( J ) . UE ( J) .RHOE(J) 
WRITE(6.2001) 

DO 211 1=1. NO 
TREF=T(I. J) *198.5 
TRRF=T(I. J)**1.5 
VIS=2.27*TRRF/TREF/10.**8. 

6= U ( I. J)*RHO ( I. J) /32.2/VIS /I 0.**6. 

UR=U( I.J) /U£( J) 

PR=P(I.J)/PH( J) 

RR=RHO(I. J)/RH0E(J) 

PTR=PTU(I.J)/PT01 
TTR=TT(I. J)/TT01 
UU=UR*UE( J) 

PO=RR*RHOE(J) 

RR0=R0/32 .2 

2S(I)=(1.-UR*RR)* (R( J)-Y(I. J) ) 

ZH(I) = (1.-UR) *UR*RR*(R< J) - Yd.J)) 

WRITE (6. 2 002) I , Y ( I. J) . YY ( I . J) . 8. EM < I. J) . UR.RR. PR. TTR. PTR 

211 CONTINUE 
SUD=0. 

SUT=0. 

DO 212 1=2. NO 
IS=I-1 
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SUO^SUO^^O.S^CYCI, JI-Y(ISf J) I^CZSC I)*2S(ISJ I 
SUr=SUT+0.5»IYII, J)-YIIS, J) >♦ CZH( I J +ZH ( ISII 
212 CONTINUE 
NO=NO*1 

TR£F=TE(J)+19a.5 
TRRF=TE(J)**1.5 
VIS=2.27»TRRF/TRE F/10.»»8« 

B=U£ CJl*RH0£CJ)/32. 2/ VIS/10. 

YZY=1. 

UR=1. 

PR=PEE(J)/PWC Jl 
RR=RHC£CJ)/RHO£CJI 
TTR=TTE(J)/TT01 
PTR=PTUEC J)/PT01 

MRITEC6, 20021 NQ* DLL ♦ YZ Yf B. EFM^ UR .RR»PR»TTRtPTR 

SUD=SUD^0.5»(OLL-Y(NOt Jll^ZSfNOI 

SUT=SUT^0.5»(DLL-y iNOfJll ♦ZHINOl 

SQ0=RCJ1*R(J1 -2.*SU0 

SQT=R<J)»R(J) -2.»SUT 

SRO=SORT(SQOI 

SRT=SORTCSQTI 

DELSTA=RC JI-SRD 

THETA=R( Jl-SRT 

WRIT£(6^2C071 OELSTAtTHETA 

WRITE (6* 2 CIO) CHILE (I) ,1-1.12) 

WRITE (6,2 000) X(J),RCJ) ,PH€J) ,OELCJ) ,CFF(J) , UE ( J) , RHOE IJ) 
WRITE(6,2006) 

WRITE(6,200i) 

CO 311 1=1, N 
Z«I)=Y(I, J) 

ZZ(I)=P(I,J) 

Jll CONTINUE 

CALL PRFL (J,TWTO,UZ,YZ,NO,UTSA) 

DO 312 1=1, NO 

T (I, J) = CTWT0<‘C1.-TWT0)^UZ4II-GAM6»EEM»EEM*UZCI)*UZ4I) / (1. 
1+GAM6^EEM^EEM) )»TT£(J) 

TE1=T£(J)/T(I,J) 

TER=SQRT(TE1) 

EMCI, J)=EEM»UZ(II »TER 
YCI,J)= 0£L(J)*YZCII 

TTCI,J) = T(I,J)^(1.*GAM6^EH(I, JI^^EM(I,J)I 
YQ=Y(I,J) 

CALL INTPCZ,ZZ,YQ,PQ,NI 
P (I,JI=PQ 

PTU(I,J)=P(I, J)*(1.*GAM6»EM(I,J)»EM(I, J))*»GAM3 
TREF=TCI. Jl+198.5 
TRRF=TCI, J)»»1.5 
VIS=2.27*TRRF/TR£F/10.»*3, 

RHOCI, J)=P(I, J)»1A4.»32.2/(1716.^T(I,J) ) 

B=UZ(I)*UH( J) ♦RHO a,J)/32.2/VlS/10.»»6. 

PR=P(I,J)/PW(J) 

RR=RHC(I, J)/RHOE( J) 

PTR=PTU(I ,J)/PT01 
TTR=TT(I, JI/TTOl 

ZS(I)=(l.-UZa)^RR)^CR(J)-Y(I,J)) 

ZH(I) = C1.-UZC I) l♦UZCI)♦RR♦(RIJ)•Y a, J) ) 

WRITE (6,2 00 2) I , Y ( I, J ) ,Y Z( 1) ,3 ,£M ( I, J ) , UZ f I) , RR, PR, TTR, PTR 


312 CONTINUE 
SUO=Q. 

SUT=0. 

DO 313 I=2»NO 
IS=I-1 

SUO=SUDi-0 J) )* CZSI IM^ZS C IS ) > 

SUr=SUT*^0 .5^(Y(I, JJ-Y (IS» Jl )» (Z HI I M^ZH (IS) I 

313 CONTINUE 
SQO=R(J>^R(JI-2**SUD 
SQT=RCJI»R(J)-2.»SUT 
SRD=SQRTCSQOI 
SRT=SQRT(SQTI 
DELSTA=R(J)-SRO 
THETA=R(J)-SRT 

WRITE! 6, 20081 DELSTA t THET A, UTSA 
8 CONTINUE 

IFIAJB.EQ.1.1 GO TO 90 
ICOD FORMATC7F10.6) 

ICIO F0RMATC12A6I 

2010 FORMAT(1H1,30X,12A6I 

2C0P FORHATC/aXdOHSTATION X^tF10.6,3H IN 1 3X » 2HR=» FIO • 3 ♦ 3H IN,3Xf 
16HPWALL=»F10.6,5H PSI A , 3X ,4H0EL= , Flu .6 , 3H IN^ 3X ,3HCF=, F10#6» 

1 2Xt3HUE=,F7.2,2X,5HRH0E=fF9^6/l 

2001 FORMATIIH * 4X » IHl » 9X , IH Yt 9X^ 2HY Y» 8 X» 6HRE < FT ) t 7X , 2HME* 9X ♦ 
14HU/UE»8X,8HRH0/RH0E»4XfifHP/PW, 8X , 6H TT/TTE 1 6 X 1 6 HPTU/PT/I 

2002 F0RHAT(2X,I3,2X* 9F12.6) 

2003 FORHAT(/aX,10HSTATION X= »F10.6»3H IN • 3X * 4HPT0=, FIO. 5, 3Xt 
14HTTO=,F10.5/) 

2004 FORMA T(/12X»lHI,12X,lHYtl2X»lHP,12X,2rtPTtl2X^lHH,12X,2HTT 

2005 FORMAT(10X,I3,2X,5F14.6l 

2306 FCRMAT(/4CXt28H WALL-WAKE VELOCITY PROFILE /I 

2007 FCRMAT(20X,7HO£LTA»=,F10.6»4H I N, , 4X ,6 HTH ET A= , FIO . 6 , 4H IN, ) 
20C8 FORMAT(20X,7HO£LTA»=,F10, 6*4H I N, , 4X .6 HTHET A= ,F 10 • 6 , 4H IN., 
14X,7HUT/U£*=,F10.6I 
END 


SUBROUTINE LWLW C J ,Z , SNN ,NO ,L, TW TOl 
DIMENSION URAT(IOO) ,Y !i00), UT (1001, PXA! 1001 
DIMENSION Z(63> 

COMMON AID,BK,CA,GAMMA,GAM1,GAM2,GAM3,GAM4,GAM5,GAH6,EEM, 

1 AN (9) ,CFFI9) ,0£L (9),RH0£ (9) ,UE C 9) , TE ( 9) , 

2 U(60,9) ,AAK,UTFT 
AN1=AN<J» 

N=AN1 

TRRF=TE(J» ♦♦l.S 

TREF=TECJM-i98.5 

VIS=2.27*TRRF/TREF/10,^»8. 

PEO£L=UE( J) »RH0E(J)*DEL ( J ) /12 .732 . 2/VIS 

SIGMA=GAM6^EEH^^2./(1.+GAM6*EEM^»2.)/TWT0 

SIGM8=l./TWT0-l. 

SA=SQRT(3IGMA> 

SAB=SIGMB^»2.+4.*SIGMA 

AS3=SGRTCSAB) 

BAS=«2,*SIGMA-SIGM8)7ASB 
DO 2 1=1, N 
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IFfZ(I)*0EL(J)l2«3v4 

2 CONTINUE 

3 NO=I 

GO TO 30 

4 NO=I-l 

30 00 1 1=1^ NO 

URATCI)=U(ItJ>/U£fJI 
Yf I)=Z(II/OEL(J) 

1 CONTINUE 
K=1 

CIAA=CA-a.614/(BK*AID) 

SIG2=2.»SIGMA/ASB 

SIG3=ASIN€BASI 

FSIG=SIG3/SA 

UT4KI=UTFT 

VWVE=C1./THT0-SIGMAI*^1.76 

AAE=REOEL*FSIG^VMV£ 

6 IF(UTCK)> 5,7,5 

7 PXA(K)=0,5 
SMM=0, 

SNN=0. 

GO TO 9 

5 PXA(K)=0.5»C1.-UHKI»((1*/BKI »ALOG C REBEL* ABS C UT (K) I ♦FS IG* VWVE ) 
l^CIAAll 

SMM=0. 

SNN=0, 

9 00 10 1=2, NO 

ACS = 1,+C0SC 3, 14*Y m ) 

YAA=SQRT{1.-YCI)**AI0I 

ALG=ALOG( Y( I) ) +2 . • ( Y A A- ALOG ( 1 . YA A 1 1 /A ID 
ACT=1 .+UT (K) ♦ALG/BK-PXA (KI*ACS 
AFF=SIN(SIG3*ACTI /SI G2 4- SI GMB/ 2* /S I GM A 

ACO=< AtG*Q*5*ACS* (ALOG( ABS(UT (KM ) ♦ ALOG ( A A£l + BK»CI A A*1 • II /BK 

SUN=(URAT (II-AFF) »FS IG* COS C SIG3*ACT) *ACD 

SMH=SMM+SUN 

AOO=ACS*UT(K)/ <2.*BK*OEL( J) )-YAA*UT<KJ /(3K*DEL <J) l-PXA «KI*Y(II* 
13*li**SINt3.14*Y (I n/DEL U) 

SUU=(URAT (I)-AFF) *CCS CSIG3*ACT) *FSIG*ADD 
SNN=SNN*SUU 

10 CONTINUE 

CFF< J)=2.*UTf Kl *FSIG*FSIG*(l./TWTO-SIGMAI *UTCKI 
IFfABSfSriM) .LE^O. COOOOOOl) GO TO 20 
IF<K,GT,i) GO TO 11 
K=K*1 

UT(K)=0.8*UT<K-1I 

SMM1=SMM 

GO TO 6 

11 K=K*1 

sl6pe=smmi/csmm-smmii 

UT(K)=UTIK-2) 4^SL0FE»(UTIK-2I-UTCK-1I ) 

S«rtl=SHH 
GO TO 6 
20 RETURN 
END 



SUBROUTINE PRFL t J ,TWTO» URT, Y, KAA, UTII 
DIMENSION URT(60) ,Y(60) 

COMMON ArO,8K*CA,GAMMA, GAHl , G AM2, GAM3t GAM4 « GAMS «GAM6«E£Mv 

1 AN (9) ,CFFC9) , DEL (9 ) t RHOE (9 ) t UE( 9) . T£ (91 t 

2 U(60*9) .AAK.UTFT 
TRRF=TECJ)»»1.5 
TREF=TE<J)+19a.5 
VIS=2.27*TRRF/TR£F/10.*^8. 

REOEL=UE( J) ♦RHOEC J>*OEL ( J ) / 12 ./ 32 • 2/ VIS 
SIGMA=GAM6*£EM**2,/a*+GAM6*E£M**2.) /TWTO 
SIGM0=1./TWTO-1. 

SA=SQRT(SIGMAI 

SAB=SIGM0»*2.+4*»SIGMA 

ASB=SQRTISAB) 

6AS=( 2.»SIGMA-SIGMB>/ASB 
CIAA=CA-0 •614/(BK»AI0I 
SIG2=2.*SIGMA/ASB 
SIG3=ASIN (BAS) 

FSIG=SIG3/SA 

UT2=CFF( J)/(2.*FSIG*FSIG*(l./TWTO-SIGMAM 
UTi=SQRr (ABS(UT2) ) 

VWVE=(l./TWTO-SIGKA) »M.76 
IF(UT1.GT*0.) GO TO 2 
PXA=0«5 
GO TO 1 

2 PXA=0.5*(l,-UTi^( (1./BK)»AL0G(RE0£L*ABS(UT1)*FSIG^VWVE) i'CIAAl) 
1 AN1=AAK+1« 

KAA=AN1 

Y(1)=0. 

URT(1)=0« 

DO 10 I=2»KAA 

AI=I-1 

Y(I) = AI/AAK 

YAA=SQRT(l.-YCI>»*AIO) 

ALG=ALOG( Y(in+2. ''(YAA-ALOG(l, + YAA))/AID 
ACS=l.+COS(3.14*Y(I)) 

ACT=1 .+UT1^ALG/BK-PXA»ACS 

URT(I)=SIN(SIG3*ACT) /SIG2+SIGMB/2./SIGMA 
10 CONTINUE 
RETURN 
END 


SUBROUTINE INTP (X ,Yt XOT tYOT.NOI 
DIMENSION X(6(3) ,y (60) ,V (60) ,VV(60 I 
DO 11 1=1, NO 
IF{X(I)-X0T) 11,12,12 

11 CONTINUE 

12 NM=I 
NU=NM+2 
NL=NH-2 

IF(NL-LT.l) NL=1 
IF(NU.GT.NO) NU=NO 
NW=NU-NL<‘1 
DO 13 J=1,NW 
L=NL*J-1 
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VIJ)=XCLI 

VViJ)=Y(L) 

13 CONTINUE 

CALL LAGRANCV^VVf XOT«yOT,NW) 

RETURN 

END 


SUBROUTINE LAGRAN (X»Y ,XOT, YOT ^Nl 
DIMENSION X«60)tYC60> 

YOT=0. 

DO 1 I-1«N 
AL=1. 

00 2 J=1»N 
lE(J.EQ.Ii GO TO 2 
AL=AL*(X0T-X« J> )/(X(I)-X(JII 
2 CONTINUE 
1 Y0T-Y0T+AL»YCI1 
RETURN 
END 


TABLE 1-B: INPUT TO PROGRAM LEAST 


0.0 

0. 

5.1 

0.4 


34.130 

1. 

1.22it05 

1.02 

1.224 

3.445 

6^243 

6.605 

8.715 

9.143 

11.215 

11.571 

12.74 

12.844 

13.065 

13.065 


TEETER 

DATA 

1. 

0 . 

1. 

539.6 

50^ 

1.9 

40^ 

4.072 

4.679 

6.922 

7.337 

9.483 

9.778 

11.851 

12.055 

12.925 

12*979 

13. 065 

13.065 


0. 

G. 

0.04 


36 « 

0.005 

5. 122 

5*542 

7.6C6 

7.959 

10.176 

10.592 

12.280 

12.462 

13.022 

13.055 

13.065 

13. 065 


0 . 


5,928 

8.350 

10.91 

12.60 

13.06 



r> r> o r> or>r>r>or»oooo 


TABLE 2-A: 


PROGRAM CONE 


PROGRAM CONE UNPUT>OUTPUT»PUNCH»TAPE 5-INPUT»TAPE 6«0UTPUT» TAPE 7 
1-PUNCH) 

INPUT FORMAT 7F10.6 
CARD(S) COLMNS 

1 1-10 GAMMA -l.A SPECIFIC HEAT RATIO OF AIR 

11-20 EMI FREE STREAM MACH NUMBER 

21-30 PHIC HALF ANGLE OF CONE TIP IN DEGREE 

31-AO XL STATION OF OUTPUT DESIRED 

THIS PROGRAM COMPUTES THE CONICAL FLOV PROPERTIES GIVEN 
THE SHOCK ANGLE AND THE SURFACE MACH NUMBER 
DIMENSION EMINF(12)»PHC(6) 

DIMENSION PH25 (12),PH50( 12 ) » PH75 ( 1 2 ) » PH100( 12) , PH125 (12 )»PH150( 12 ) 
DIMENSION EM25a2)>EM50(12)/EM75(12)»EM100(12) »EM125(12)»EM150(12) 
DIMENSION PC(6)»0(6) 

DIMENSION SAVE1(300)> SAVE 2(300 ),S A VE3(300),SAVE<» (300 )» SAVE? (300) 
1»SAVE6(300)»SAVE7(300) 

DIMENSION SAVEBOOO) 

TT-2000, 

R-53.3 

THE FOLLOWING CAROS READ IN THE DATA AS TAKEN FROM THE 
CONICAL FLOW TABLES BY SIMS 

DATAEMINF/1.?»1.75»2.C»2.5>3.0»3.5»A.C»A.5,5.0/6.0»7.0»8.0/ 

DATA PHC /.0<t 3633 ».0e726b» .13090, .17A 53 » .21 817» .26180/ 

DATA PH25/.72980,.60832,.5237C,.A1169, .3A011,. 29017, .25329,. 22A95, 
1. 2C252,. 1693 7, .1A623, .12931/ 

DATA PH 50/. 73706,. 60953,. 5 2525,. 41 425, .34910,. 29592, .26104, 

1.2348 4, .214 58, .18565, .16626, .15262/ 

DATA PH75/. 73 490,. 61445,. 5 31 40,. 42358,. 3 5701,. 31239,. 28078, 

1.25747. . 23973, .21481, .19646, .18712/ 

DATA PHIOO/. 74466, .62564, .54465, .44136, .37900, .3379C, .30918, 

1.28622.. 27242. .25049. . 23633. .22664/ 

DATA PH12 5/. 76 158, .64 399,. 5651 8, .4 6631, .40760, .36939, . 34296, 
1.32363, .30953, .2 8990, .2 7740, .26895/ 

DATA PH150/.78592,.669C0, .59192, .49663, . 44085, .40494,. 3803 2, 

1.36266. . 34955. .3 3174, .32052, .31301/ 

DATA EM25/ 1.48669, 1.73403, 1.98086, 2. 47295, 2. 96 276, 3. 45009, 
13.93477,4.41662,4.89553,5.84407,6.77969,7.70188/ 

DATA EM 50/ 1.4 5793, 1.7004 6, 1.94 152, 2.41940,2.89138,3.35716, 
13.81645,4.26906,4.71481,5.58525,6.42692,7.23908/ 

DATA EM75/ 1.41 973,1.65680, 1.89123,2.35286,2.80480,3.24675, 
13.67839,4.09946,4. 50971, 5.29688,6.03849,6.73371/ 

DATA EMIOO/ 1.374 86, 1.6064 5, 1.8 3404, 2. 27 6 88, 2. 7101, 3. 12751, 
13.53057,3.91891,4.29218,4.99278,5.63160,6.21039/ 

DATA EM125/ 1.32490, 1.551 33, 1.772 19, 2. 200 16, 2. 6 10 39, 3. 00266, 
13.376145,3.73078,4.06627,4.68017,5.22045,5.69227/ 
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DATA EM150/1.2707<r,l.A9255»1.7C688#2.1179<*^2. 50675/ 2. 87312, 
13.21670,3.53743,3.83559,^.36681,^.81780,5.19815/ 

1 READ (5.1010J GAMMA, EMI, PHIC, XL 
WRITE(6,?000) GAMMA, EMI, PHIC, XL 
PHIC -PHI C/57. 2957795 
KKK-1 

GAMMl-GAMMA-1.0 

GAMP1«GAMMA+1.0 

THE NEXT PART OF THE PROGRAM INTERPOLATES TO FIND THE 
SHOCK ANGLE AND THE SURFACE MACH NUMBER 


DMON-0.0 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
AMDN-0.0 
CALL TAINT 
DMON-0.0 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
CALL TAINT 
BMON-0.0 
CALL TAINT 


(EMINF,PH25,EM1,PC{1),12,2,NER,DM0N) 
(EMINF,PH50,EM1,PC(2),12,2,NER,OMON> 
(EMINF,PH75,EMl,PC(3), 12, 2,NEP,0M0N) 
(EMINF,PH100, EMI, PC (4 ) , 12 ,2, NEP, DMON ) 
t EMINF,PH125,EM1,PC(5), 12,2,NER,DMOM) 
(EMINF,PH150, EM1,PC(6),12,2,NER,DM0N) 

(PHC,PC,PHIC,PHIW,6,2,NER,AM0N> 

(EMINF,EM25,EM1,0(1),12,2,NER,DM0N) 
(EMINF,EM50,EM1,Q(2),12,2,NER,0M0N» 
(EMINF, EM7 5,EM1,0(3),12,2,NER,DM0N) 
(EMTNF,EM1C0, EM 1, Q{ 4) , 1 2, 2, NER, DMON ) 
(EMINF,EM125,EM1,C(5J,12,2,NER,0M0N) 

( EMINF,EM150, EM1,Q(6),12,2,NER,DM0N) 

( PHC,Q,PHIC,EMC,6,2,NER,BM0N» 


calcutate the maximum velocity 


20 VM-S0RT(2.0/GAMHl)*S0RT(GAMhA*R*TT*32.2» 


initialize the velocity in the angular direction 

3C PSIO-0.0 


CALCULATE THE VELOCITY ON THE CONE 


40 VROU"( EHC**2)* (VM**2)*(GAMMl/2.0) 
VROD«1.0+(EMC**2)*(GAMM1/2.0» 

VR0»SCRT(VR0U/VP0D) 

PHICC-PHIC 

C CALCULATE THE ENTROPY FUNCTION 

D1»AL0G( t <7.0)*(EM1**2)*(SIN (PHIW)**2»-l.C>/6. C) 
D2-( <6.0)4 (EM14*2) ♦(SlN{PHIW)**2l) 

D3»( ( EM1442)*(SIN(PHIW)4*2)+5.0) 
0S-(Dl-(1.4)*AL0G<D2/03) )/0.4 
DP*EXP(-0S) 

SAVE1(KKK)«XL 

SAVE2(KKK )»XL*S IN (PHIC) /COS ( PHIC) 
SAVE3(KKK)«PHIC*57. 2957795 
SAVE4(KKK)»EMC 
SAVE5(KKK)«DP 

S A VE6(KKK)-PHIC4 57. 2957795 
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Ill III ■! 


1 1 ■■ll■ll■llll■l■| ■ I III I I II I 


ll■■lll II 


mini mil I i 


III ■■■■iiiii I III iiimi II 


EMX-EHC 

A1-1.+0.5*6AMM1*EM1*»2. 

A11-AL0G(A1» 

B1*A11*GAMMA/GAMM1 

Bll-EXP(Bl) 

A2«1.+0.5*GAMH1+EHX**2. 

A21»AL0G( A2) 

B2*A21*GAMMA/GAMH1 

B21*EXP{B2) 

SAV£7(KKK)«DP*P11/B21 

KKK-2 

C SET THE INCREMENTAL VALUE OF THE ANGLE TO BY USED 
50 DELPHI* (PHIW-PHIC) /117.0 
WRITE I6>2050) 

60 DO 200 K*l»117 
STEP-K 

X-STEP*OELPHI 
PHl-PHICC+X 
FUZZ-STEP-1.0 
IF(FUZZ) 68»68»70 

C CALCULATE THE VELOCITY IN THE AN6UALR DIRECTION 
68 P3*(GAMM1)*{VR0)*( ( VM** 2 ) - ( V R0**2 )) 

D«-(GAMMl/2.0) ♦( ( VM*+2)-(VR0**2») 

PSIX-PSI0+(DELPHI+P3)/D 

GO TO eo 

70 P1»-(GAMM1/2.0)*(PSI0**3)*(C0S(PHIC) /SIN(PHIC) ) 

P 2 »(GAMMl/ 2 . 0 )*( (VH**2)-( VR0**2n*(PSI0)*(CQS(PHIC)/SIN(PHICn 
P3*(GAMK1)*( VRO)*l ( VM**2 ) - ( VR0*+2 ) ) 

PA»-(GAMMA)*( VRO)*( PSI0**2 ) 

D=(GAMPl/2.0)*(PSTO++2)-(GAMMl/2.0)*((VM**2)-( VR0**2» ) 
PSIX*PSI0+((DELPHI)i'(PH-P2 + P3 + PAn/D 
C CALCULATE THE VELOCITY IN THE RADIAL DIRECTION 
80 VRX*PSIO*DELPHI+VRO 
C CALCULATE THE FLOW ANGLE 
100 THETAX*PHI+ATAN(PSIX/VRX» 

110 VX*SORT (VRX++2+PSI X**2) 

120 AX*SORT( (GAMMl/2.0)*( VM**2-VX**2) ) 

C CALCULZTE THE MACH NUMBER 
130 EMX-VX/AX 

EMlU*(EMX**?)+(GAMPl/2.0) 

EM10«1.0+( EHX**?)*(GAMMl/2,0> 

EM1ST»SOPT(EM1U/EM10 ) 

140 YL«XL*{SINIPHn/COS(PHin 
BURP-PHI+57. 2957795 
WRITE (6»2100) VRX,PSIX.BURP»AX 
WRITE (6,2110) XL, YL,THETAX,EM1ST,DP 
SAVE1(KKK)»XI 
SAV€2(KKK)»YL 

SAVE 3 ( KKK ) .THETA X+57. 295 7795 
SAVE4(KKK ) -EMX 
SAVE5IKKK) -OP 
SAVE6(KKK)-BURP 
A1=1.+0.5*GAMM1*EM1**2. 

A11-AL0G( Al) 

B1-A11*GAMMA/GAHM1 

Ell-EXP(Bl) 

A2-1.+C.5*GAMM1*FMX**2. 


I 
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A21«AL06( A2) 

B2>A21*GAMHA/6AHH1 
B21«EXP(B2 ) 

SAVE7(KKK)-0P*R11/B21 

SAVEe(KKK)»1.02-SAVE2(KKK) 

PSIO-PSIX 
VRO-VRX 
PHIC-PHI 
KKK-KKK+1 
200 EMC-EMX 

WRITE(6»3000» 

WRITE (6, 3001) ( S AVE1(KKK),SAVE2(KKK)»SAVE3(KKK)»SAVE4(KKK)» SAVES (K 
IKK), SAVE6(KKK)»SAVE7(KKK),KKK«1,11 8»2) 

WRITE! 7,3002) ( SA VE6 (KKK ) , S AVE3 (KKK ) ,SA V EA( KKK ) , S A VE7( KKK) , KKK*1, 
1118,2) 

WRITE(7,3003) (SAVE1(KKK),SAVE2(KKK),SAVE3(KKK),SAVE4(KKK), 
1SAVE5(KKK), KKK-1, 118,2 ) 

GO TO 1 

3000 FORMAT! 1H1,12X,1HX,1AX,1HY,IAX,3HDEL,13X,1HM,15X,8HPT/PTTNF,10X, 
15HTHETA,12X,AHP/P1) 

3001 FORMAT! 2X,7F15. 6) 

3002 FORMAT! AFIO. 6) 

3003 F0RMAT!5F10.6) 

1000 FORMAT !?F10.3) 

1010 FORMAT !AF10,6) 

2000 FORMAT I1H1,7HGAMMA «, F 1 5 .6, lOX, AHMC «, F 15. 6, 1 CX , 6HPHIC •» 
1F15.6,10X,3HX •,F15.6) 

2050 FORMAT ! IHl, ASX, 7HRESULTS ) 

2100 FORMAT !1H0,AHVR » , F 1 5 . 6, 1 OX, 5HPS I •,F15.6,10X, 

16HPHI »,F15.6,10X,3HA »,F15.6) 

2110 FORMAT !1HC,3HX • , FIO . 5, 8X, 3HY •, FI 0 . 5, 8X , 8HTHET AX -,F10.5,8X,* 
16HMSTR -,F10.5,8X,10HPT/PTINF «,F10.5) 

END 


SUBROUTINE TAINT!XTAB,FTAB,X,FX,N,K,NER, MON) 
DIMENSION XTAB !l),FTAB! 1),T!10),C!10) 

CPSOAOO TAINT SUBROUTINE- IN FORTRAN II, 

IF !N - K) 1,1,2 

1 NEP-2 
RETURN 

2 IF !K-9) 3,3,1 

3 IF ! MON) A, A, 5 

5 IF ! MON-2) 6, 7, A 
A J«0 
NMl-N-1 
DO 8 I»1,NM1 

IF !XTAB!I )-XTA5! I+l) ) 9,11,10 
11 NER-3 
RETURN 
9 J-J-1 
60 TO 8 
10 J«J+1 
8 CONTINUE 
MON-1 

IF !J) 12,6,6 
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12 M0N«2 

7 DO 13 I»1>N 

IF {X-XTAB(D) 1A,14,13 

1^ j-r 

GO TO 18 

13 CONTINUE 
GO TO 15 

6 DO 16 T*1>N 

IF (X-XTAB(D) 16^17#17 

17 J-I 

GO TO 18 
16 CONTINUE 
15 J-N 

18 J«J-(K+l)/2 

IF (J) 19fl9^20 

19 J*1 

20 M*J+K 

IF (M-N) 21^21^22 

22 J-J-1 

GO TO 20 

21 KP1«K^1 
JSAVE-J 

26 DO 23 L-1#KP1 
C (L )»X-XTAB( J) 

T(U*FTAB( J) 

23 J-J+1 

DO Zk J-1>K 
I-J + 1 

25 T(n*{C(J)*T(I)-Cm + T(J))/(C(J)-C(D) 
I-I + l 

IF (I--KP1) 25f25 
2^ CONTINUE 
FX-T(KPl) 

NER-1 

RETURN 

END 


TABLE 2-B; INPUT TO PROGRAM CONE 


l.A 2.82 10. 1.835 
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TABLE 3-A 


PROGRAM ANAL 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM ANAL ( INPUT , OUTPUT » TAPE 5= INPUT , TAP£6=0UTPUTI 
ANA 

TWO-DIMENSIONAL AND AXIALLY SYMMETRIC SHOCK HAVE-BOUNO ARY LAYER 
INTERACTION 


PROGRAM INPUT 

INPUT FORMAT 7F10.6 EXCEPT CAROS 1 AND 2 
CARD(S> 

1 AND 2 TITLES, COLUMNS 1-72, HOLLERITH 

3 COL.l-iO EMEl UPSTREAM MACH NO. 

11-20 UTUESl UPSTREAM DIMENSIONLESS FRICTION VELOCITY 
=UE/UE^ 

21-30 REOELl UPSTREAM REYNOLDS NO. (DELTA) 

31-40 AINPT l.C INPUT THETAl 

2.0 INPUT ALPHAl 
41-50 RECFAC RECOVERY FACTOR 

51-60 AXI 0.0 OR 1.0 2-D CASE 

2.0 AXI CASE CONE FLOW OR CENTERBODY 
61-70 SHORT 0.0 LONG OUTPUT 

1.0 SHORT OUTPUT 

4 COL. 1-10 BLEEDl SUCTION RATE 

11-20 BLEED2 ENTRAINMENT RATE 

21-30 AIB 0.0 NO BLEED 

1.0 POROUS SUCTION 

2. C SLOT SUCTION 

3. C SCOOP SUCTION 

31-40 8K CONSTANT IN WALL-WAKE PROFILESCUSUALLY 0.4) 

41-50 C CONSTANT IN WALL-WAKE PROFILES (USUALLY 5.11 

51-60 STCOEF 0. IF NO SHEAR FORCE IS DESIRED 

0.5 IF AVERAGE SHEAR STRESS IS USED 
61-70 AIOO PARAMETER IN ( 1 .- ( Y/ DELT A) AIDD) ( = 1. IS 
RECOMMENDED AT THIS TIME) 

5 COL. 1-10 THID THETAl WHEN AINPT=1. 0 

ALPID ALPHl WHEN AINPT=2.0 
11-20 TH1Z2 THETAl WHEN AINPT=2.0 
FOLLOWING FOR AXIALLY SYMMETRIC CASES ONLY (AXI=2.0I 

6 COL. 1-10 OELIR DELTAl/R 

11-20 AI£2 1.0 B.L. EDGE STREAMLINE DATA INPUT 

2.C CONICAL FLOW DATA INPUT 

3.0 INTERNAL FLOW FIELD INPUT DATA 

4. C CONSTANT PRESSURE BOUNDARY 

21-30 FLDIR3 FLOW DIRECTION ONSTREAM OF AXISYM. INTERACTNS 

7 COL. 1-10 AKK NO. POINTS OF INPUT DATA 

FOLLOWING FOR B.L. EDGE STREAMLINE I NPUT < AIE2=1. 01 

8 COL. 1-10 XXO X/XO 

11-20 RRO R/RO 

21-30 EME2S MACH NUMBER 

31-4C PE2P1 P/Pl 


I 
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C 41-5D ALPHA FLOW DIRECTION 

C CA1?D 6 IS REPEATED AKK TIMES IN ORDER FROM SHOCK TO CONE 
C FOLLOWING FOR CONICAL FLOW I NPUT <A IE2=2. 01 
C 8 COL. 1-10 PHI CONICAL RAY ANGLE 

C 11-20 ALPHA FLOW DIRECTION 

C 21-30 EME2S MACH NUMBER 

C 31-40 PE2P1 P/Pl 

C CARD 8 IS REPEATED AKK TIMES IN ORDER FROM SHOCK TO CONE 
C FOR SUCTION WITH NO SHOCK USE INPUT CAROS 1-6 WITH 
C AINPT=2.0, ALP10=0.0 
C 

c 

DIMENSION TTRATC200) ,PTRAT«200I ,PTRNS( 200 ) ♦URATI 200 ) ^EM (2001 
DIMENSION YY(20f»* UTUE3(100)t EIIOO), WR(200J» PHIR(200) 
DIMENSION REOELSdOOl 


COMMON 

TTRAT 

f 

PTRAT 

t 

PTRNS 

• 

URAT 

t 



1 

£M 

♦ 

WR 

f 

PHIR 

t 





2 

SHORT 


BLEED 1 

t 

8LEED2 

» 

BLEED 

• 



3 

AIB 

» 

TWITE 


GAHMA 

9 

THIO 

f 




BLMN 


BLMON 


8LMI 

9 

BLMOI 

• 



5 

OSO 


□□SO 


BK 

9 

EMEl 

9 

SHFAC 

9 

6 

SIGl 

« 

SIGMAl 

f 

SIGSl 

9 

FSIGl 

9 

VWVEl 

9 

7 

CFWWl 

t 

C 

f 

STCOEF 







COMMON /UWU/AIDD 
COMMON /WUW/AINPT^THIZZ 
C 

C READ DATA CAROS 1 AND 2 
C 

1 CALL PRIN (ILlNEtll 
C 

C READ DATA CAROS 3 AND 4 
C 

R£AD(5»10C0I EMEl»UTU£Sl»REDELltAINPTtRECFAC,AXI,SHORT 

READ (5. 13 00) BLEEDl,BLEE02»AIB,BKtCtSTC0EF^AI00 

C=C-0.614/C3K^AIOOJ 

IKPT=AINPI 

BLE£D=BLEE01«-BLEE02 

GAMMA=1.4 

G1=EME1*EME1 

THTTE=R£CFAC+ (1.- FECFAC)/ (1.+ IGAMMA-1. ) »Gl/2.^ 

SIG1=0.2»£ME1»EME1 

SIGMA1=SIG1/(1.+SIG1) 

SIGS1=SQRT(SIGMA1) 

FSIG1=(ASIN(SIGS1) l/SIGSl 
VHV£1=(TWTTE*(1.+SIG1II ♦♦l./e 
FUTUE1=UTU£S1*FSIG1 

CFWW1=2.»FUTU£1^FUTUE1/ (T HTTE»( 1. ♦$! G1 ) I 
TAU1P1=0. /♦CFWWl^EM£i*EM£l 

PXl=.5»(l.-UTUESl»((l./eK)^AL0G (REOELl* ABS C UTUESl) 
l»FSIGl/tfWV£l)+C)> 

C 

C VARIABLES WHICH MUST BE INITIALIZED 
C 

IZ=0 

YS=0. 

WRITE(6,1024) 

GO TO (1Q«20I »INPT 
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INPUT PARAMETER IS THETAl 
READ DATA CARD 5 

10 REAOC5»1000) THiO 

WRITE (6»10i0) UTUeSltR£0£LifPXl,THlD,EMEl»RECFAC*T«TTE 

ILINE=ILIKE+3 

TH1=TH10». 0174532 9 

Z=SINCTH1)**2. 

ANUM=2.»1G1»Z-1.) 

DN0M=(2.frGl^(GAMMA+l.-2.»Zn»TANCTHll 

ALP1=ATAN<ANUM/0N0MI 

ALP10=ALPl/#ai745329 

FIND PRESSURE RATIO ACROSS SHOCK AND DOWNSTREAM MACH NO* 

CALL SHOCK(2t ALPl«fTHl*EMEltEME2*P2Pll 
GO TO 30 

INPUT PARAMETER IS ALPHAl 
READ DATA CARO 5 

20 READC5. 10001 ALPlOfTHlZZ 

WRITE C 6, 1011) UTUES1*RE0EL1»PX1,ALP10,EME1.RECFAC*TWTTE 

ILINE=ILINE+3 

ALP1=ALP1D^. 01745329 

IFCALPIO.EQ.O*) GO TO 30 

FIND SHOCK WAVE ANGLE, PRESSURE RATIO ACROSS SHOCK AND DOWNSTREAM 
MACH NO. 

CALL SH0CK<1,ALP1,TH1,EME1,EME2,P2P1I 

TH1D=TH1/. 01745329 

ROOT=SIN(THl) 

ERROR EXIT IF SINCTHETA) GREATER THAN 1.0 

IFCROOT.GT.l. JGO TO 1 

EQUATE ALPHA! AND ALPHA2 

30 ALP2=ALP1 
ALP20=ALP10 
L=1 

WRIT£(6, 10251 L,8LEED1 
L=2 

WRIT£(6,lC25i L,BLE£02 

FIND SHOCK WAVE ANGLE, PRESSURE RATIO ACROSS SHOCK AND DOWNSTREAM 
MACH NO. 

IFCALPIO.EQ.O.) GO TO 310 

CALL SHOCK (1,ALP2,TH2,EM£2,£HE3,P3P2I 

ROOT=SIN(TH2l 

ERROR EXIT IF SINCTHETA) GREATER THAN 1.0 
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IF(R00T,GT*1.IG0 TO 1 
TH2D=TH2/.01745329 

310 IFCALPID-GT.O.) GO TO 311 
C 

C PAPAMETERS FOR NO SHOCK 

C 

£H£2=EME1 

EME3=EME1 

P2P1=1. 

P3P2=1, 

F3P1=1. 

GO TO 312 

311 P3P1=P3P2»P2P1 

WRITE C6» 1012) P2Pl»P3P2,P3Pl,EM£l,ErtE2,EM£3 
ILINE=ILINE^►6 

WRITE (6,1013) ALPlOfTHlO, ALP20,TH2D 
ILINE = ILINE+6 
C 

C DETERMINE UPSTREAM BOUNDARY LAYER PROFILE PROPERTIES 
C 

312 CALL PRFL(UTUES1,PX1,£ME1,1.0,2,YY) 

IF(SHORT,EQ.O.) CALL PRIN ( ILl NE ,3 ) 

11=1 

WRirE(6,1014l II 
ILINE=lLINE+4 

IFCSHCRT.EQ.l.) GO TO 4C00 
WRirE(6,1028) 

DO 4C 1=1,101 

WRITE (6, 1015) I,y Y(I> ,EM(I) , URAT ( I ) , PTRAT ( I ) , PTRNS ( I) , TTRATC I ) 

ILINE=ILINE+1 

CALL PRIN(ILINE,2) 

IF(ILINE-2) 31,31,40 
31 WRITE(6,1014) II 
WRITE(6,1028) 

ILINE=lLINE+4 
40 CONTINUE 
4000 OUH=0. 

C 

C DETERMINE UPSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 
C 

CALL FLUX(101,YY,1,EME1,DUM) 

HI1=(SHFAC-SIG1)/ (THTTE»(1.+SIG1) ) 

RETH1=REO£L1*OOSO 

PRANDL=*72 

PR13=PRANDL*».ll./3.) 

T3ARTC=.5»TWTTE^^* 22»PR1 3+ ( . 5- •22»PR13) / (1 . +SIG1 ) 
CFLT1=«246»EXP(-1 .561»HI1I » (RETH1»* (-.268)1/ 

1( (TBARTO* (l.^-SIGl )) **0. 79 631 
WRITE (6,1 016) DSD,0DSD,BLMN,BLM0N,SHFAC,HI1,CFWW1,CFLT1 
D5D1=0SD 
DDS01=00SD 
ILINE =ILINE*6 
C 

C CONSTANTS FOR MASS AND MOMENTUM BALANCES 
C 

Hl = l*+ (GAMMA- 1.1/2. ♦EME3^EM£3 
H2=SQRT(H1) 
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H3=1.*CGA«MA-1#)/2.»G1 

Hh^SQRT(H3I 

H5=1,-DSO-OOSO 

IFCAIB.EQ.l.) GO TO 410 

IFC3CEED1.EQ«0.I GO TO 410 

WS=ABS (BLEEOli ♦BLHN 

CALL lNTRP(WR»YYf WS,YSf 1011 

CALL INTRP(YY»PHIR,YS,PHIS»101I 

WRIT£«6,1C27) WS,YS,PHIS 


AIB=1 POROUS PLATE SUCTION 
AIB=2 SLOT SUCTION 
AIB=3 SCOOP SUCTION 


ITERATION TO FIND SOLUTION TO CONTINUITY ANO MOMENTUH EQUATIONS 


410 IFCAIB.LE.l*! BLMOMl^O* 

IF(AIQ-GE*2.) BLMCH1=-GAMMA»ENE1»EME1»PHIS 
KB=1 

WRITE (6» 10261 K0*BLMOM1 

BLHOM2=3LEE02 ♦ I 1 • -0 SO ) »G AMMA*EN£l»E ME2^C0S ( ALPl) 

l*H4/SQRTf 1.^IGAMMA-'1*)/2.»EM£2*EME2> 

KB:=2 

WRITE(6f 1026) KB» BLH0M2 

BLM0MF=BLM0rtH>BLM0M2 

KB=3 

WRIT£C6*1026I KBtBLHOMF 

FOR APPROXIMATE SOLUTION ASSUME THAT TAU/P1=2 . ♦STCOEF^TAUl/Pl 
TAUP1=2.»STC0EF»TAU1P1 

H6=P2P1-1.-GAMMA*G1^H5-BLM0MF+TAUP1/TAN(ALP1I 
IF(AIB.EQ.3.> H6=H6- (P2P1-1.) ♦YS 

IF CSHORT.EQ.a .1 W FITE ( 6 , 13 2 2) HI * H2 1 H3 ,H4 » H5 1 H6 »T AUPl 

Al=P3Pl^£ME3/£MEl*H2/H4/4i,-OSO)/(l.«-BLE£OI 

A2=(P2P1-P3P1)/H6 

A3=-GAMMA*P3P1*£ME3*EME3/H6 

A4=TAUP1/TAN(ALP1)/H6 

IFCSHORT.EQ.O.) WRITE <6 , 1 017) Al,A2tA3tA4 
IFISHORT.EQ.O.) CALL PR IN 4 ILI NE »3 ) 

IF(SHORT.EQ.O.I WRITE <6 ♦ 1 01 dl 
ILINE=ILINE<^2 

APPROXIMATE 2-D SOLUTION CA SSUME RE0EL3=R£DEL1 AND TAU=TAU1) 


M=1 

R£D£L3(1)=REDEL1 

SIG3=.2»£M£3^EME3 

SIGMAS^SIGS/d.^-SIGS) 

SIGS3=SQRT(SIGMA3) 

FSIG3=4ASIN43IGS3))/SIGS3 
VHV£3= (TWITE* 41.+SIG3 ))♦»!. 76 
UTU£341)=1.5*UTUES1 
IFIAIB.GE.l.l UTU£3(1I=3,*UTUES1 

PX3=.5*(1#-UTU£341)* ( ll./8K)*ALOG (RE0EL3 4M)*ABS (UTUE3 ( 1) ) 
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1*FSIG3/VMVE3) ^-Cll 
J=1 

GO TO 43 
41 J=J*1 

IFCJ^EQ.aC) GO TO 50 
UTU£3<JI=UTUe3(J-ll-«01 

PX3=.5'»^(1.-UTU£3(J)» (Cl./8K)*AL0GfR£D£L3«M)»ABS(UTUE3« Jll 
l*FSIG3/\/WVE3) +Cn 

43 CALL PRFL(UTUE3(Jl5PX3,£ME3,P3Plt ItVyi 
CALL FLUX(10i,YY,l»EME3»DUMI 
E{J) = A2<-A3-AH‘f A1-A3)*OSD-A3^DOSO<-A4 
IF(SHORT. EQ.l.) GO TO 440 
WRIT£l6tl019l J*UTU£3<JI,ECJ) 

ILlNE = ILlNEi-3 
CALL PRIN (ILINEt2} 

IF(ILINE-2) 430f430f440 
430 WRITE(6tl018l 
ILINE=ILINE+2 
440 T£ST==ABS(EC Jl) 

IF(TEST.LE*0,000C1I GO TO 60 
IF(£(J).LT,0*l GO TO 45 
IF<J,LT.9> GO TO 41 

45 SLOPE=(E( J-1) -E(J>J/(UTUE3(J-H-UTUE3f Jll 
UTUE3 <J>1 )=UTU£3i Ji-£ (Ji/SLOPE 
47 J=J^1 

PX3=.5*(1.-UTUE3( Jl»f (1 ./BK)*AL0G(R£0EL3(M)*ABS(UTU£3I JM 
1^FSIG3/VWVE3) «-C) > 

IF(J*EQ.8C) GO TO 50 
GO TO 43 
50 WRITEC6,1C20I 
12=1 
GO TO 80 


SOLUTION 08TAINE0 (APPROXI MATE! 


DETERMINE DOWNSTREAM BOUNDARY LAYER PROFILE PROPERTIES 

6C CALL FRFL (UTU£3CJ),PX3t £M£3,P3Pit2^YYI 
IF(SHORT*EQ*0.) CALL PRIN ( ILI NE t 3> 

11=3 

WRITE(6tl014l II 

WRITE<6,1033I UTUE3(JJ, REDEL3 C Ml ,PX3 
lLlNE=lLlNE+6 

IF(SHORT.EQ.l.i GO TO 7000 
WRITE<6, 10281 
DO. 70 I=ltl01 

WRITE(6,1G15I I,YY(II ,EM(II »URAT< II , PTRAT ( I) , PTRNS ( 1} , TTRATI I) 
ILINE=ILINE*-1 
CALL FRIN <ILIN£,2I 
IF(ILIN£-2) 61,61,70 
61 WRITEI6,1C14) II 
WRITE(6, 10281 
ILINE=ILIN£+4 
70 CONTINUE 
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DETERMINE DOWNSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 

7000 CALL FLUXClOlf YYtl»EME3»0UMJ 
FUrU£3=UTUE3(JI^FSIG3 

CFWW3=2**FUTUE3»FUTU£3/ <TWTTE*(1. ^'SIGSII 

FH1=1.+SIG1 

FM3=1.+SIG3 

HI3=(SHFAC-SIG3»/(TWTTE»FM3I 

RETH3=R£0EL3<M)*0DSD 

TBART C=.5^TWTTEfr.22*PR13M#5-.22*PR13l /I1.+SIG3I 
CFLT3=.246»EXP(-1.561»HI3J*«RETH3»*C-.268) 17 
l((T8ART0*(i.*^SIG3n ♦♦0.79631 
WRITE (6,10161 OSDtDOSO, BLMN , BLMON , SHFAC tHI3 , CFWW3 , CFLT3 
D13=A1^(1.-0S01 
031=1. /D13 
DS31=0SD^031/0S01 
DOS 3l=D0SD»D31/ODSOl 
IF(ALP10.EQ.O.) GO TO 79 

INTERACTION LENGTH 

EL01= (l.-C3i) /TAN (ALPl) 

IFCAIB.EQ. 3.1 £L01=(1.-031-YSI/TAN<ALPH 
79 IF(ALP10.EQ.Q.) ELOl=0. 

WRITE (6,1021) D31,US31,D0S31,EL01 

TWO-DIMENSIONAL EXACT SOLUTION (ONLY FOR 2-0 CASES) 

IFCM.GT.l) GO TO 80 
IF(AXI.EQ.2.) GO TO 80 
IF(031.LT.O.) WRITE (6,10341 
IF(031.LT.C.) GO TO 1 
CALL PRIN(ILINE,3) 

WRITE(6,1C29) 

IF (SHORT. EQ.l.) GO TO 500 
WRITE (6, 1030) 

WRITE (6,10311 M,RE0EL3(M) 
lLINE=lLlNE+4 
500 M=MH 

IF (H.EO.ll) GO TO 530 

RE0HL3(M)=R£DEL1^D31»P3P1»CEM£3/EME1)^(FM3/FM1) *^1.26 
FUTUE3=uruE3( JI^FSIG3 

CFWW3=2.»FUTUE3^FUTUE3/ (T WTTE»( 1. +SIG3) ) 

TAU3P1=3. 7^CFWW3^EME3^EME3»P3P1 
TAUP1=STC0EF^(.TAU1P1«*TAU3P1) 

IFCSHORT.EQ.l.) GO TO 650 
WRIT£(6,1031) M,REDEL3(M) 
lLlNE=ILlNE+2 
CALL PRIN(ILINE,2) 

IFIILINE-2) 600,600,650 
600 WRIT£(6,1030) 
lLlNE=ILlNE+2 

650 TESTR= A3S ( ( RE0EL3 CM-1) -RE0EL3 CM ) I /RE0EL3 (M-1 J I 
IF(TESTR.LE..C01) GO TO 60 
UTUE3 (1)=1.5^UTUES1 
IFCAIB.GE.l.) UTUE3(1)=3.^UTUES1 

PX3=,5»(1.-UTUE3(1)^( (1./8K)»AL0G(RE0EL31M)»A8S(UTU£3C1)) 
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1*FSIG3/VMVE3>»C)I 

J=1 

GO TO 903 

5C1 J=J+1 

IFCJ.EQ.OOl GO TO 90 
UTUE3(J) = UTUe3(J-l)-,.01 

PX3=.5»il*-UTUE3IJ>*l a./BK)»AL0G<R£0EL3IM>*AaSCUTU£3IJII 
1»FSIG3/VMV£3I+CII 

503 CALL PRFL(UTUE3tJ)»PX3»£M£3,P3Pl,l,YY> 

CALL FLUX (101, YY»l,EME3tDUMI 

H6=P2Pl-l.-GArtMA»Gl^H5-BLMOHF+TAUPl/TAN(ALPll 
IF(AIB.EQ.3.) H6=H6- (P2P1-1.I*YS 

IF(SHORT.£Q*0.) W RITE ( 6 , 102 2) HI, H2, H3 , H4 , H5,H6 ,T AUPl 

A2=(P2P1-P3P1)/H6 

A3=-GAMMA*P3P1*£ME3*EM£3/H6 

A4 IS THE SHEAR TERM IN THE MOMENTUM EQUATION 
A4= TA UPl/TANC ALPD/H6 

IFCSHORT.EQ.O.) W.RI TE (6 , 10 17) A1,A2,A3,A4 
E( JI=:A2^A3-AH-(A1-A3)*DSD-A3*D0S0^A4 
TEST=ABS(E(J)I 

IF(T£ST.L£.0.D0001) GO TO 520 
IF(J.LT.3) GO TO 501 

5 05 SLOPE=(£l J-1) -E(J))/(UTUE3(J-1) -UTUE3(J)I 
UTUE3(J4*1) = UTUE3(J)-£(J) /SLOPE 

PX3=.5*(1.-UTUE3( J) ♦((!./ BK)*ALOGIREDEL3( ri)*ABS(UTUE3 (J) I 
1^FSIG3/VWVE3I*C)> 

IF(J*EQ.8C) GO TO 50 
GO TO 503 

520 D13=A1»(1.-0S0I 
D31=l,/D13 
GO TO 50 0 

530 WRITE(6,1C32I 
IZ=1 

GO TO 60 

60 IF(AXI.LT.2.) GO TO 1 

IFCALPl.GT.O.I CALL AXI S CELDl ,UTUES1,UTU£3I J) , 
1RE0EL1,PX1,ALP1,EME1,I2) 

IFCALPl.EG.O.) CALL ANOSCUTUESl ,UTUE3< J) , REOELl,PXl, 
lD31,EMEl,IZi 
GO TO 1 

ICOO FORriAT(7F10*6) 

1010 FORMATCIHO, 9X,2lH0PTI0N 1 INPUT THETAl/1 CX ,9HUTUES1 = ,F10.6, 
15X,24HREYN0L0S NO. (DELTAl) = , E14 .6 ,5 X, 6HPX1 = ,Flfl.6/5X, 
25X,9HTHETA1 = ,F10.6,BH D EGREES ,5 X, 5HM1 = ,F10.6/ 
310X,16HRECOVERY FACTOR = , FIO . 6 ,5X , 9HTW/TTE = ,F10.6) 

1011 FCRMATCIHC, 9X,2lHOP7ION 2 INPUT ALPHA 1/1 0 X , 9HU TUES 1 = -,£10.6, 
15X,24HREYN0L0S NO. (DELTAl) = , £14. 6,5 X,6HPX1 = ,F10.6/5X, 
25X,9HALPHA1 = ,F1C.6,8H DEGREES ,5 X,5HM1 = ,F10.6/ 
310X,18HRECOV£RY FACTOR = , F 10 .6 , 5X, 9HT W/T TE = ,F10.6) 

1012 FORMAT(1HO/13X,3CHBOUNDARY LAYER EDGE CONDITIONS// 

11QX,8HP2/P1 = ,F10.6,5X,8HP3/P2 = , F 10 . 6, 5X , 8HP3/P1 = ,F10.6/ 
210X,5HM1 = ,F10.6,8X,5HM2 = , FI 0. 6, 8 X, 5HM3 = ,F10.6) 

1013 FCRMAT(1HO/1CX,32HSHOCK ANO FLOW DEFLECTION ANGLES//10X, 
115HINCIDENT SHOCK ,5X,8HALPHA = ,F10.6,8H DEGREES, 5X, 
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28HTHETA = ,F10,6, 6H OEGREES/lQX »15HREFLECTE0 SH0CK,5X^ fiHAtPHA * « 
3F10.6»8H 0EGR£ES*5X,8HTHETA = ,F10.6,8H DEGREESI 

1014 FORHATUHC, 9X, 35H80UN0 AR Y LAYER PROFILE DATA STATION, 131 

1015 FCRHATCIH ,6X, 13, 6F14.6 > 

1016 FORMATClHO/10X,14HOELSrAR/QEL = , F14 . 6 ,4X , lOHHOM/DEL = , 
1F14.6/1CX,28HN0N-0IMENSI0NAL MASS FLUX = ,F10,6/ 
2lOX,32HMON-DIMENSIONAL MOMENTUM FLUX = ,F1Q.6/1QX, 

315HSHAPE FACTOR = , FlO. 6, 5X,3CH INCOMPRESSIBLE SHAPE FACTOR = t 

4 F10.6/10X,16HCF(WALL-WAKEI = ,F10 .6,5X,22HCFCLUOW£IG TILLMAN) « , 
5FiO,6l 

1017 FORMAT C1H0,10X,5HA1 = ,Fia • 6,5X ,5 HA2 = ,F10 ,6 ,5X,5HA3 = ,F10,6, 

1 5X.5HA4 = ,F10.6) 

1018 FORMAT<1HO/1CX,32HRESULTS OF ITERATIONS FOR UTUE^S//) 

1019 FORMATC1HC,10X,4H J = ,I 3, 5X ,9HUTUE^3 = ,F1 0 .6 ,5 X,4H£ = ,F10.6I 

1020 FORMAT ClHO, 6X,25HNO CONVERGENCE FOR UTUE»3) 

1021 FCRMAT<1HC/10X,12HOEL3/DEL1 = ,F10.6/10X, 

120HDELSTAR3/OELSTAR1 ~ ,F10. 6/lOX ,12HMOM3/MOM1 = ,F1C.6//10X# 
29HL/3EL1 = ,F10.6) 

1022 FORMAT(1HO/10X,5HH1 = , FlO . 6, 2X ,5HH2 = ,F10 .6 ,2X,5HH3 = ,F10,6, 
12X,‘5HH4 = ,F10.6,2X,5HH5 = ,F10 • 6 ,2X ,5HH6 = ,FlO-6/ 

2 1QX,9HTAU/Pt = ,F1Q.6) 

1024 FCRMAT(lHO/iOX,54HTMO-DIMENSIONAL APPROXIMATE SOLUTION (REDEL3 = R 
lEOELin 

1025 FORHATCIHO, 9X , 3H (M ) BLEED , II, 11 H/ (M> B. L* = ,F10,6, 

120H (BLEED POSITIVE IN)) 

1026 FORKATdHC, 9 X , 40HMOMENTUM FLUX ASSOCIATED WITH B.L. BLEED, II, 

13H = ,F14.6,14H (POSITIVE IN) I 

1027 FORMATC1HO/10X,5HWS = , FlO • 6, 5X ,5HYS = ,F10.6, 5X,7HPHIS s , FlO. 61 

1028 FCRHATdHC, 

1 9X,1HI,8X, 1HY,14X,1HM, 12X,4HU/U£, 9X,6HPT/PT1, 6X , ICHPT 4 NSI /PTl, 

2 €X,6HTT/TT1//I 

1029 FORMAT (1HO,7X,30HTWO-DIMENSIONAL EXACT SOLUTION) 

1030 FORMAT(1HO/10X,31HRESULTS OF ITERATION FOR RE0EL3I 

1031 FORMAT(1HO,10X,4HM = , I 3, 5X , 22HRE YNOLD S NO. (DELTA) = ,£14.61 

1032 FCRMAT(lHOt6X,25HNO CONVERGENCE FOR REOEL31 

1033 FCRMAT(1H0,10X,9HUTUE^3 = , FlO . 6, 5X, 9H REDEL3 = ,E14.6,5X, 

16HPX3 = ,F10.b) 

1034 FORMAT ( 1 HO , 6X , 28HOEL3/OEL1 NEG. (NO SOLUTION)) 

END 


SUBROUTINE AXIS(£L01,UTU£S1 ,UTUE3 ,RED£L 1,PX1, ALP1,EME, IZI 
C AXIALLY SYMMETRIC FLOW ANALYSIS 

DIMENSION TTRAT(200l,PTRAT(20a) ,PTRNS ( 200 ) ,URAT (200 ) ,EM 1200) 
DIMENSION £14100), E2 4100 ) ,EL01 AdOQ ) , UTUE3A 4100) , YY (200) 
DIMENSION XXO( 50), RRO ( 50),£ME2S( 5u ) , PE2P1C 50 
DIMENSION PHI4 50, ALP HA 4 50), RR ( 5G), DEL30R4 50) 


DIMENSION EI2( 

50 

), £L02( 

501 ,XR4 

50) 

,WR(200) , 

PHIR4200) 

COMMON 

TTRAT 

f 

PTRAT 

, 

PIRNS 


URAT , 


1 

EH 

t 

WR 

t 

PHIR 




2 

SHORT 

9 

BLEEDl 

, 

BLEE02 


BLEED • 


3 

Aia 

9 

THTTE 

, 

GAMMA 


THID , 


4 

BLMN 

9 

BLN ON 

, 

BLMI 


BLMOI , 


5 

DSD 

, 

OOSO 

♦ 

9K 


EMEl , 

SHF AC 

6 

SIGl 

f 

SIGMAl 

♦ 

SIGSl 


FSIGl , 

VWVEl 

7 

CFWWl 

t 

C 

• 

STCOEF 
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COMMON /HUW/AINPTtTHlZZ 


READ DATA CARO 6 

PEAD(5,10C0) DEU1R,AIE2,FLDIR3 
CALL FRFL(UTUESl»PXl,£H£l,1.0t2fVY) 

CALL FLUX (101, YYt 2»EHE1»0£L1RI 

SN1=C1^-0SD)/D0SD 

SN2=SNl-3*0 

SN3=AL0G(SN2> 

SN4=-0.6169*SN3 

SN5=£XP(SN4) 

F£N=0*0299^SN5 

F£NT=F£N»ELD1/BLMN 

FEHT=F£NT*BLMI/BLKOI 

PG0=1.-DEL1R 

ROOl=RO0*ROO-2.*BLMI»FENT 

R0R1=SQRT(R001) 

I£2=AI£2 

ALP10= ALP 1/. 01745329 
CALL FRIN(ILINE*3) 

WRir£(6,lCll) 

GO TO (10,20,35), IE2 

BOUNDARY LAYER EDGE STREAMLINE DATA ARE INPUT 
READ DATA CARO 7 

10 PEAD(5, 10001 AKK 
KK=AKK 

DO 11 1=1, KK 
READ DATA CARDS 8 

READ (5, 100 0) XXO(I) ,RRO (I),EME2S(II ,PE 2P1 < I ) , ALPHA ( II 

11 CONTINUE 
WRITE( 6,1025) 

ILINE = ILlNE^-5 
DO 12 1=1, KK 

WRITE(6,1C26) I , XXO( I) , RRC ( I) ,EME2S1I) ,PE2P1( I) , ALPHA! I) 

ILINE=ILINE+1 

CALL PRIN CILIN£,2) 

IF(ILINE-2) 13,13,12 
13 NRITE(6,1025I 
ILINE=ILINE«-5 

12 CONTINUE 
GO TO 25 

PROGRAM CALCULATES CONICAL FLOW STREAMLINE FROM INPUT DATA FROM 

CONICAL FLOW PROGRAM 

ARRAYS IN ORDER SHOCK TO CONE 

20 R£AO(5,1000) AKK 
KK=AKK 

DO 21 1=1, KK 
PHI IS IN DEGREES 
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ALPHA IS IN DEGREES 

REAO<5»1300> PHKII, ALPHAdI i EHE2S<II» PE2PU19 

21 CONTINUE 

IF<SHORT,£Q,0.> WRITE (6 ♦ 10271 

ILlN£=lLlN£+5 

DO 822 I=lfKK 

PHI ( I > =PH I ( IJ ♦ . 91745329 

ALPHA < II = ALPHA < II *•01745329 

PHI IS IN RADIANS 
ALPHA IS IN RADIANS 

IF(SHORT.£Q.l,l GO TO 822 

WRITE 16* 10281 I*PHI C 1 1 * ALPHA! II *EME2SC 1 1 »PE2P1I IJ 

ILINE=ILINE+1 

CALL PRINIILINE»21 

IF(ILIN£-2J 621^821*822 

821 HRITEC6* 10271 
ILIN£=ILINE«-5 

822 CONTINUE 

CONICAL FLOW STREAMLINE CALCULATION 

TERM1=9. 

TERM2=0. 

XXO!ll=l* 

RR0fl)=l. 

TP1=TANCPHI(1)I 
DO 22 K=2»KK 

PHIBAR=(PHI(K-1I+PHI(K) 1/2. 

ALPBAR=(ALPHA(K-ll+ALPHA!KII/2. 

TPa=TAN(PHlBARI 
TAB=TAN( ALPBAR) 

TERM1 = T£RM1 + I1,*TPB*TP8J* (PHUKI-PHI (K-ll ) / IT AB-TPB) 
XX0(KI=£XP(T£RH1) 

TERH2 = TERM2dAB*(XX0 ( K) -XXO IK-1 ) I 
RR0(K)=1. +TERM2/TP1 
IFCXXOlO .GT.20.1 GO TO 922 

22 CONTINUE 

922 IF(SHORT*EQ.O.I WRITE 16 • 1 025) 

ILINE=ILINE«*5 
OO 24 I=1^KK 

ALPHA C II = ALPH A ( II / . 0 1745329 
ALPHA IS IN DEGREES 

IFCSHORT.EQ.1.1 GO TO 24 

WRITE(6,1026) I* XX0III» RROIlIt EME2S C II , PE2P1 1 1 ) , ALPHA I II 

1LINE=ILIN£+1 

CALL FRIN (ILINE*2) 

IF(ILIN£-2I 23,23,24 

23 WRITEI6,1025) 

ILINE=ILIN£^5 

24 CONTINUE 

CONVERSION OF STREAMLINE COORDINATES TO RS/R ANO L/DELl 
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25 R0R=1.-0EL1R 
ROR=ROR1 

TH 11= TH13^0. 01745329 

IF(AINPT.EQ.2, ) TH11=TH1ZZ»0. 31745329 

TP1=TAN(TH11I 

X0R=1.»R0R/TP1 

DO 26 K=1»KK 

RR(K)=RR0 (K)»R0R 

XR(K)=XX0 <K)*R0R/TP1 

EL02(K)=(XR(K)-X0RI/0EL1R 

DEL30R(K)=1.-RR«K) 

26 CONTINUE 

IFiSHCRT.EQ.l.) GO TO 2630 
WRITE(6,1029) 

ILINE=ILIN£^5 
DO 28 1=1 ,KK 

WRITE (6,1033) I , XXO (I ) , RRC( I) ,RR(II ,EL02(I) ,0EL30R ( I) , XRf I) 

ILINE=ILIN£+1 

CALL PRIN (ILINE,2I 

IFCILINE-2) 27,27,28 

27 WRITE(6, 10291 
ILINE=ILINE+5 

28 CONTINUE 

CALCULATION OF FORCE ALONG BOUNDARY OF CONTROL VOLUME ADJOINING 
REGION 2 

2800 EI2C1)=0. 

DO 280 K=2,KK 
RRBAR=(RR (K-1) ^RRIKD/E. 

PE2Plfi = CP£2Pl(K-lH*P£2Pl(K) )/2, 
EI2(K)=EI2(<-1)-PE2P1B*RRBAR*(RR(K-1)-RRIKII 
280 CONTINUE 

IF(SH0RT.EQ.1« ) GO TO 3000 
WRITE(6,1030) 

ILINE=ILINE+5 
DO 30 1=1, KK 

WRITE (6, 1031) I,XXO(I),RRC(I) ,RR( I) , ELD2( I ) ,0EL30R (I) , £NE2S( I ) 
1,PE2P1(I), EI2CI) ,ALPHA(I) 

ILINE=ILIN£<‘1 
CALL PRIN(ILINE,2) 

IF(ILIN£-2) 29,29,30 

29 WRITE(6,1C30) 

ILINE=ILINEfr5 

30 CONTINUE 
3000 GO TO 37 

CONSTANT PRESSURE BOUNDARY 

FIND SHOCK WAVE ANGLE, PRESSURE RATIO ACROSS SHOCK AND DOWNSTREAM 
MACH NO. 

35 CALL SHOCKd, ALP1,TH1,£ME1,EM£2,P2P1I 
ALP1D= ALP 1/. 01745 329 
ALP2=ALP1 
ALP20=ALP10 
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C FIND SHOCK HAVE AHGLE« PRESSURE RAriO ACROSS SHOOT AN 0 OOWNSTREAM 
C MACH NO. 

C 

CALL SHOCK II, ALP2»TH2,£ME2,EMEJ»F3P2B 
TH2D=TH2/. 01745329 
P3P1=P3P2*P2P1 
WRITE(6tlQ34l 

WRITE 16, 10121 P2P1,P3P2,P3P1,EM£1,EME2#B«3 
ILINE=ILINE<*6 

WRITE (6, 1013) ALPlOfTHlDf ALP2d,TH20 
ILINE=ILINE+6 

DETERMINE UPSTREAM BOUNDARY LAYER PROFILE PROPERTIES 

37 CALL PRFL(UTU£S1,PX1,EME1,1.0,2,YY) 

IFISHORT.EQ.O.) CALL PRIN(ILINE.3) 

II==1 

WRITE(6,1G14) II 
ILIN£=ILINE+4 

IFISHORT.EQ.l.l GO TO 4000 
HRITE(6, 10391 
DO 40 1=1,101 

WRITE 16, 1015) I,YY(I) ,EM(I) ,URATC I) ,PTRAT 1 1) , PTRNS ( I) , TTRATi I) 
ILINE=ILINE+1 
CALL FRINIILINEtE) 

IF(ItINE-2) 38,36,40 

38 WRITEC6,1014) II 
ILINE=ILINE+4 
WRITE(6,1C39) 

40 CONTINUE 

DETERMINE UPSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 

400D CALL FLUX(101,YY,2,EH£1,DEL1R) 

HI1=I SHFAC-SIGl)/ (THTTE»I1.+SIG1) ) 

RETH1=REOEL1*ODSO 
PRANDL=.72 
PR13=PRANDL*» (1./3.) 

T3ARTC=.5^TWTTE<'.22^PR13^'I.5-.22*PR13)/(1,*SIG1) 
CFLT1=.246*£XP(-1.5&1*HI1) *{RETHl*^l-.268) )/ 

1( CTBARTO^ (1.+SIG1))*»C.7963) 

WRITE (6,1016) OSD,ODSD,BLMN,BLMON,SHFAC,HI1,CFHW1,CFLT1 

ILINE =ILIN£+6 

DS01=OSO 

DOSD1=ODSO 

WRITE (6,1017) D£L1R,BLHI,8LM0I 
ILINE=ILINE+3 
1 = 1 

UTU£3ACI)=UTU£3 

FM1=BLMI 

FM01=BLMOI 

FM01=FM01»(1.+FEMTI 

FM1=FM1»(1.+FENT) 

J=1 

ELDIAC J)=EU01 

R3=BLEeD*2,»FMl 

L=1 
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WRITE (6»1C36> LtBLEEDl 
L=2 

WRITE(&tlC36l t,6L££D2 
WRITEC6»1035I BLEE0fR3 
IFCIZ.EQ.il RETURN 
IFCEL01 .lt. 0.1 RETURN 
YRS=0. 

IFCAIB.EQ.l.) GO TO 500 
IFC3LEED1.EQ.0.I GO TO 500 
WSA=ABSC 3LEED1) *2.»FM1 
CALL INTRPCWR.YY, WSA,YYS»101I 
CALL iNTRPCYY^PHIR^YYSf PHISA,101I 
YRS=YYS»0EL1R 

WRITEC6tl0371 WSA tYRS »PHISA 
R3R=1. 

IFCAie.EQ.3.1 R3R=1.-YRS 

AIB=1 POROUS PLATE SUCTION 
AIB=2 SLOT SUCTION 
AIB=3 SCOOP SUCTION 

500 IFCAIB.LE.l.) BLMCM1=C. 

IFCAIB.LE.3.) R3R=sl. 

IFCAI8.GE.2.I BLMCMl^-PHISA 
K8=l 

HRITEC6»1038i KB*BLMOMl 
IFCSHORT.EQ.1.1 GO TO 50 
CALL PRINCILIN£,3I 
WRITE(6tl018l 
ILINE=ILI^E♦^ 

50 GO TO C51*51»51t55lt IE2 

ITERATION PROCEDURE FOR SOLUTION OF CONTINUITY EQUATION 

51 CALL INTRP{£LD2,EME2S,£LD1A(J1,EME2#KKI 
CALL INTRP{EL02tPE2Pl*ELDlACJ),P2PltKK) 

CALL lNTRPC£L02.ALPHAtEL01ACJ) ,ALP2D,KKI 
ALP2=ALP20».017<»5329 

ALLOW FOR NON-ZERO FLOW ANGLE DOWNSTREAM* FL0IR3 

ALP2=ALP2-FL0IR7». 01745329 
ALP20=ALP2/. 01745329 

DETERMINE SHOCK HAVE ANGLE* PRESSURE RATIO ACROSS SHOCK AND 
OOWNSTREAM MACH NO. 

CALL SH0CKC1*ALP2*TH2*EME2*EME3*P3P2I 
P3P1=P3?2+P2P1 
TH2D=TH2/. 01745329 
IFCSHORT.EQ.l.) GO TO 55 

WRITE (6* 1012) P2P1*P3P2*P3P1*EME1,EHE2*EME3 
ILINE=ILINE+6 

WRITE (6, 1013) ALP ID * THIO* ALP20.TH20 
ILINE=ILINE+6 
55 Rl=2. 

GAHl=IGAMMA-l.)/2. 
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R1A=1 • f'GAMl*EME3^EM£3 
RlB=l.*GAMi*EMEl*EMEl 
R2=-2,*P3P1^EME3*SQRT<R1A/R1BI/EME1 
IFIAIB.EQ.3.) R2=R2*R3R»R3R 
R3=BLEED»2*»FH1 

IFCIE2.EQ.3) D31P=1.-ELD1AC J> *TANCALP1I 

IFIIE2.EQ.3.AN0.AIB.EQ.5. ) 031P=1 . -YYS-ELDl A < J) *T AN (ALPII 
IFCIE2 .lt. 3> CALL INTRP ( ELD2, DEL30R, EL DIA C J) » DEL3RP,KKI 
IFCIE2.LT.3.AN0.AIB.EQ.3.) 0£L3RP=DEL3RP-YRS 
IF(IE2.LT.3» 031P=DEL3RP/CEL1R 

REDEL3=RED£L1»D31P*P3P1»CENE3/EME1)*(R1A/R13I»»1.26 

SIG3=GAM1*EME3^EME3 

SIGMA3=SIG3/R1A 

SIGS3=SQRT(SIGMA3I 

FSIG3=CASIN(SIGS3))/SIGS3 

VWVE3=ITWTTE*R1A>*»1.76 

IFCREOEL3.LE.O.) RETURN 

PX3=.5*I1.-UTUE3A Cl ) * ( 1 1. /BKI *ALOG (RE0EL3* A9S <UTUE3A<I)) 
l^FSIG3/YWV£3)+Cn 

60 CALL PRFL(UTU£3A(II ,PX3,EME3»P3Plfl.YY> 

GO TO (61, 61.61, 62» , IE2 

61 call INTRP(£LD2.EI2,ELD1ACJ>.EI2PR,KK> 

CALL INTRP(EL02.DEL30R,ELD1A( Jl .DEL3R.KKI 
IFCSHORT.EQ.l.) GO TO 63 

WRITE (6.1032) ELOIA C J) . EI2PR. DEL3R 

ILINE=ILINEH 

CALL FRIN(ILINE.2) 

IF(ILINE-2) 610,610.620 
610 WRITE(6. 10161 
ILlNE=ILINE+2 
620 GO TO 63 

62 D31=1.-EL01A( JI^TAN(ALPl) 

DEL3R=D31»DEL1R 

63 C0N=D£L3R 
D3PR3=(DEL3R-YRS) /R3R 
IF(AIB.EQ.3.I OEL3R=D3PR3 
IFCSHORT.EQ.O.) WRITE (6 .1040) R3R.D3PR3 
lLlNE=ILIN£+2 

CALL FLUX (1Q1.YY,2.EN£3.0EL3R) 

0EL3R=C0N 

FM3=BLMI 

FMO3=0LMOI 

E1CJ)=R1»FM1+R2*FM3+R3 
IFCSHORT.EQ.l.) GO TO 65 

WRITE(6. 1C19) Rl. R2, R3. FMl. FHOl ,F M3. FM03 
WRITE(6.1020I I.UTUE3A(I) .E2(I) . J ,EL01A ( Jl .El CJ) 
ILINE=ILINE+6 
CALL PRIN (ILINE.2) 

IFCILlNE-2) 64.64,65 

64 WRITEC6. 10181 
ILINE=ILIKE4-2 

65 TEST=ABS(E1CJ)> 

IFCTEST.LE. 0.00001) GO TO 80 
IFCJ.GE.2I GO TO 70 

J=2 

ELDlACJ)=l.i»£L01ACl) 

GO TO 50 



UOO ouo oou uoo 


70 SLOPEl=(El{J-l)-EliJI )/ (£L01ACJ-H-ELD1A<JH 
ELOIA (JH) = ELD1A< JI-EKJI/SLOPEl 
0£FG=ELDlA(J*i) 

EFG=EU01A (1»*1.6 

IF(OEFG.GT,EFG» ELOIA (J»1I=1.1*EL01A<J» 

J=JH 

IF(J.EQ.51) GO TO 71 
IFCELOIAIJ)) 72,72,50 

STOP CASE NO CONVERGENCE ON J 

71 WRITE(6,1021» 

RETURN 

STOP CASE L/DEi.1 LESS THAN 0. 

72 MRITE(&,1CL2) J,EL01A<J>,I 
RETURN 

ITERATION PROCEDURE FOR SOLUTION OF MOMENTUM EQUATION 

80 Slsl. 

S2*-l. 

S3=-2. 

S4= 2.*GAMMA*eNEl*EMEl 
S5=2.*GAMMA*EME3*£ME3*P3P1*C-1.» 

S6=GAMMA*EME1*EME l*BLMOMl 
S7=-l. 

TAU1P1=0. 7»CFMW1*ENE1*EME1 
FUTUE3=UTUE3A(IJ*FSIG3 

CFMW3=2.*FUTUE3*FLTUE3/ tTHTTE*<l.*SIG3ll 

TAU3P1=0.7*CFHH3*£ME3*EM£3*P3P1 

TAUP1=STCCEF*(TAU1P1+TAU3P1» 

58 IS THE SHEAR TERM IN THE MOMENTUM EQUATION 

S8=2.*TAUP1*0EL1R*EL01A (Jl 
R1C=1.*GAM1*EME2*EME2 

59 IS THE ENTRAINMENT TERM IN THE MOMENTUM EQUATION 

S9=2, ♦BLEED 2*FM1*C0S«ALP1>*EME1*E ME2*GAMMA*SQRT iRlB/RlCI 

T1=2.*0£U1R-DEL1R*0EL1R 

T2=2.*O£L3R-OEL3R-»0£L3R 

IF <IE2.£G,3.) T1=2.*D£L1R+0EL1R*0EL1R 

IF (IE2.EG.3.I T2=2.*D£L3R*DEL3R»0EL3R 

T3=P2Pl»(Tl-T2>/2. 

T4=0. 

IF(AIB.EQ.3.» T2= (2. *D3PR3-03PR3*03PR3) *R3R*R3R 
IFCAia.EQ. J.» T4=2.*yRS-yRS*TRS 
IFCAIB.EQ.3.) S5=S5*R3R*R3R 
IF(IE2.LT.4I T3=EI2PR 

E2tIl = Sl*Il+S2*T2-*P3Pl«-S3*T3+Sif»FM01«-S5*FM03+S6*S7*T4-S3*S9 
IF (SHORT. EQ.l.) GO TO 82 
WRir£(6,1022> S4,S5.SS,T1,T2,T3,T4 
HRITE(6.1Q2a) I,UTUE3A(I) ,E2(I) ,J,EL01A(J) ,E1(J) 
IL1NE=ILINE*6 
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CALL PRIN(ILIN£»2) 

IFIILINE-21 81t«l,82 

81 WRITEI6»1018) 

ILIN£=ILIN£*2 

82 TEST =A83CE2Cin 
IFfTEST,LE,O.OOC01) GO TO 110 
IF<I.GE.2) GO TO 90 

1=2 

UTUE3ACII=0.9^UTUE3 

IF«EME1*LT.2.5.AN0.ALP10*LT.3.5I UTUE3»e»^l*UTUE3 
GO TO 100 

90 SL0PE2=CE2CI-1I-E2(II )/(UTUE3AII-l>-tfniE5MD I 
UTU£3A(I+1)-UTUE3AIII-E2CII/SL0PE2 
1 = 1+1 

IFCI.EQ.51) WRIT£(6^1023) 
ifci.eq.51) return 

100 ELD1AC1)=ELD1AIJ) 

J=l 

GO TO 50 

SOLUTION TO CONTINUITY AND MOMENTUM EQUATIONS MAS BEEN OETERHINEO 

110 UTUE3=UTU£3ACIJ 
EL01F=EL01ACJ) 

DETERMINE DOWNSTREAM BOUNDARY LAYER PROFILE PROPERTIES 

CALL FRFL (UTUE3 ,P X3 t FME3 f P3P1 »2 ^YYI 
IFCSHORT.EQ.O.) CALL PRIN C I LINE »3) 

IFCSHORT.EQ.a.) WRITE C6.1018) 

IFCALPIO.EQ.O.) GO TO 1100 

WRITE C6t 1012) P2Pl*P3P2,P2Pl»EMEl»EM£2tEME3 
WRITE (6, 1013) ALPlOtTHlO, ALP2D,TH2D 
1100 IFCSHORT.EQ.O.) CALL PR IN C I LINE ,3 1 
11=3 

HRITEC6tl014) II 

WRITE C6» 1043) UTUE3t RE0EL3» PX3 

ILlNE=ILlN£+6 

IFCSHORT.EQ.l.) GO TO 1200 
WRIT£C6*1039) 

DO 120 1=1^101 

WRITE C6» 1015) I*YYCI) »EMCI) ,URATCI) , PTRAT C I) , PTRNS C I) ♦ TTRA T Cl) 

ILINE=XLINE+1 

CALL PRINCILINE.2I 

IFCILINE-2) 111*111»120 

111 WRITEC6,1014) II 
ILINE=ILIhE+4 
WRITEC6,1039) 

120 CONTINUE 

DETERMINE DOWNSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 

1200 IFCAIB.EQ.3.) 0£L3R=03PR3 

CALL FLUX ClOlf YYt 2,£N£3,0EL3RI 
FUTUE3=UTUE3+FSIG3 

CFMW3=2.*FUTUE3+FUTU£3/ C T WTTE* C 1. +SI G3) ) 
MI3=CSHFAC-SIG3)/CTWTTE+C1.+SIG3) ) 



RETH3=RE0EL3»00S0 

T8ART0=.5»TWTTE^'.22»PR13+(.5-.22*PR13i/ll*i‘SIG3l 
CFLT3=.246»EXP(-1.561»HI3l»<RETH3»»f-.26dl 1/ 

If CTBARTO* Cl. + SIG31) ♦♦0.79631 
WRITE C6t 10161 0SD»DDSD,BLMN*3LM0N,SHFAC»HI3*CFWH3fCFLT3 
IF(AIB.EQ.3.) DEL3R=D3PR3^R3R 
WRITEf6tl0171 D£L3R»BLMI,BLH0I 
ILINE=ILIK£^3 
D31=OEL3R/OEL1R 
DS31=0S0^031/0SD1 
DOS 31=D0SD»O31/0OS01 
WRITE (6,1 024) 031 , DS31, D0S31,EL01F 
IF(Aie.EQ.3.) WRIT£(6fl041) YRS ,D3PR3»R3R 
RETURN 

1000 FORMAT (7F10. 6) 

ICll FCRHAT(lHC/10X»35HSOLUTlON FOR AXIALLY SYMMETRIC FLOW) 

1012 FORMAT (IH /lO X, 30H80UNOAR Y LAYER EDGE CONDITIONS// 

113X*8HP2/P1 = ,F1C.6,5X,8HP3/P2 = , FIO . 5X , 8 HP3/P1 = ,F10.6/ 
210X,5HM1 = ,F10.6,8X,5HM2 = , Fl Q* 6, 8 X, 5HM 3 = ,F10#6) 

1013 FORMAT(lHO/iaX,32HSHOCK AND FLOW DEFLECTION ANGLES//10X. 
115HINCIDENT SHOCK ,5X,8HALPHA = ,F10.6,8H DEGREES, 5X» 

26HTHETA = ,F13.6, 8H OEGREES/ICX , 1 5HREFLECTED SHOCK , 5X, 8HALPHA = ^ 
3F10.6,8H DEGREES, 5X,8HTH£TA = ,F10.6,8H DEGREES) 

1014 FCRMATflHO, 9X , 35 H80UNU AR Y LAVER PROFILE DATA STATION, 13) 

1015 FORMATdH , 6X , 13, 6F14 .6 ) 

1016 FCRMAT(lH0/10X,14HOELSTAR/DEL = , F14. 6, 4X, lOHMOM/OEL = , 
1F14.6/10X,28HNON-DIMENSIONAL MASS FLUX = ,F10.6/ 
2lOX,32HNON-OIMENSIONAL MOMENTUM FLUX = ,F10.6/10X, 

315HSHAPE FACTOR = ,F1G. 6, 5X,30HXNCOMPRE SS I BLE SHAPE FACTOR = , 

4 Fia. 6/10 X,16HCF( WALL-WAKE) = , FIO . 6,5 X, 22HCF (LUDWEIG TILLMAN) = • 
5F10.6) 

1017 FORMATdH , 9X,8H0EL/R = ,F10 .6/1 OX, 21HB. L. MASS INTEGRAL = , 
1F10.6/10X,25H3.L. MOMENTUM INTEGRAL = ,F1Q.6/) 

1018 FCRMATdHC, 9X , 42 HRESULTS OF ITERATION FOR UTU£»3 AND L/OELl) 

1019 FORMATUHO, 9X,5HR1 = , FIO . 6, 6X ,5HR2 = ,F10 .6 , 5X ,5HR3 = ,F10.6/ 
110X,6HFM1 = ,F10. 6,5X,7HFM01 = ,F10 . 6/ 1 OX , 6HFM3 = ,Fi0.6,5X, 
27HFM03 = ,F10.6) 

1020 FCRMATflH , 9X,4HI = , 13, 3X, 9HU TUE^3 = , FIO . 6 ,5X,5HE2 = ,F10.6/ 
110X,3HJ =,I4,5X,7HL/01 = , F 10 . 6 , 5X , 5HE1 = ,F1D.6) 

1021 FCRMATflHO, 6X,19HN0 CONVERGENCE ON 3) 

1022 FORMATflHC, 9X,5HS4 = , FIO . 6, 5X ,5 HS5 = ,Flfl .6, 5X, 5HS8 = ,F10#6/ 
110X,5HT1 = ,F10.6,5X,5HT2 = • FIO. 6,5X,5HT 3 = , FIO . 6,5X ,5HT4 = , 

2 FIO. 6) 

1023 FORHATflHO, 6X,19HN0 CONVERGENCE ON I) 

1024 FCRMAT(iHQ/lCX,12HD£L3/0ELl = ,F10.6/1CX, 

120HO£LSTAR3/DcLSTAR1 = ,F10 .6/1 OX ,12HM0M3/M0M1 = ,F1C.6//10X, 
29HL/DEL1 = ,F10.6) 

1025 FORMAT (1HO/10X,35HBOUNOARY LAYER EDGE STREAMLINE DATA// 

1 9X,1HK, 8X,4HX/XO,10X,4HR/RO,11X,3HM2E,10X,5HP£2P1, 

2 9X,5HALPHA//) 

1026 FCRMATCIH , 6X , 13, 5F14 .6 ) 

1027 FORMAT(1HO/10X,34HCONICAL FLOW STREAMLINE INPUT DATA// 

1 9X,1HK,9X,3HPHI,10X,5HALPHA, 9X,3HM2E, 8 X , 5HPE2P1//) 

1028 FORMATCIH , 6X , 13, 4F14 .6 ) 

1029 FORMATflHO/10X,25HCONVERTED STREAMLINE DATA// 

1 9X,lHKtaX,4HX/XO,lCX,4HR/RO,iOX,4HRS/R, 9X,6HL/0EL1, 8X,6HDEL3/R 
2,10X,3HX/ R//) 
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1030 FCRMAT(1H0/10X»31HSUMMARY OF EDGE STREAMLINE DATA// 

1 9X»lHKf4Xt4HX/XO»10X<4HR/RO, 8X,4HRS/R, 7Xf 
26HL/DEU1, 7Xt6H0EL3/Rt 8Xt3HM2£, 7Xt5HPE2Pl» 7X, 

35H(IIE2» 7X,5HALPHA//> 

1031 FORMATC IH •EX, 13 , 9F12 • 61 

1032 FCRMATCIH 1 10 X 1 9HL/0EL1 = » F14. 6f 5X» BH C I) E2 = tF14.6^5Xt 
16H0EL3R = tF14«6) 

1033 FCRMATCIH , 6X » 13* 6F14.6 I 

1034 FCRMATClHa/10X*38HCONSTANT PRESSURE BOUNDARY IN REGION 21 

1035 FORMATCIH * 9X * 19H CM ) BLEED/ (MI B .L . = *F10,6/ 
11QX*19H(MIBLE£D/(M)0UCT .F10,6J 

1036 FORMATIIH « 9X * 8H (Ml BLEED *11* IIH/ (Ml B.L. = fF10.6» 

120H (BLEED POSITIVE INI I 

1037 FORMATClH0/10X*6HWSA = ,F10 *6 *5X* 6HYRS = ^ FIO .6 * 5X » 8HPHISA = * 
1F10.6I 

1038 FCRMATCIHO* 9X » 40HMOI1ENTUM FLUX ASSOCIATED WITH B.L* BLEED*I1* 

13H = *F14.6,14H (POSITIVE INJJ 

1039 FORMATdHOt 

1 9X*lHl*aX,lHY*14X*lHM* 12X.4HU/UE* 9X*6HPT/PT1* 6X * ICHPTC NSI /PTl # 

2 ex*6HTT/TTl//l 

1040 FCRMATdHfl* 13X*7HR3/R = , FIO . 6, 5X ,12H0EL3PR/R3 =- *F10.6I 

1041 FCRMATdHQ* 8X,28HSC00P CASEtSCOOP HEIGHT/R = ♦FlO.btSX* 
112HDEL3PR/R3 = *F 1 0 . 6 *5 X* 7HR3/R = ,F10.6I 

1042 FORMATdHO/ 6X*17HST0P CASE L/D £L 1 (* I3,4HI = , F 10 . 6 *5X ,4H I = *131 

1043 FORMATdHO, 10X*9HUTUE*3 = * FI 0* 6* 5X * 9HREDEL 3 = *E14«6*5X* 

16HPX3 = *F10.6I 

END 


SUBROUTINE AN03 (U TUESl * UTUE3, REOELl * PXl * 03 12D , £M£ * IZI 
C AXIALLY SYMMETRIC FLOW ANALYSIS BLEED WITH NO SHOCK WAVE 

DIMENSION TTRAT(2 00) ,PTRAT(20C) *P TRNS ( 20C I * URAT (20 0 I * EM ( 2001 
DIMENSION El(lOO), E2 ( 1 00 I * 0EL30R ( 10 C ) * UTU£3A(100I* YY(200I 

DIMENSION WR(200I* PHIR(200I 


COMMON 

TTRAT 

• 

PTRAT 

9 

PTRNS 

9 

URAT 

9 



1 

EM 

* 

WR 

f 

PHIR 

9 





2 

SHORT 

* 

BLEED 1 

» 

BLEED 2 

9 

BLEED 

9 



3 

AIB 

• 

TWITE 

• 

GAMMA 

9 

THIO 

9 



4 

BLMN 


3LM0N 

9 

BLMI 

9 

BLMOI 

9 



5 

OSD 

'f 

OOSD 

9 

8K 

9 

EMEl 

9 

SHFAC 

9 

6 

SIGl 

V 

SIGMAl 

9 

SIGSl 

9 

FSIGl 

• 

VWVEl 

9 

7 

CFWHl 

* 

C 

9 

STCOEF 







READ DATA CARO 6 


REAO(5*1000I OELIR 
CALL PRIN(ILINE*3) 

WRltE(6*1011l 

DETERMINE UPSTREAM BOUNDARY LAYER PROFILE PROPERTIES 
37 CALL PRFL(UTU£S1,PX1*EH£1*1*0*2,YY> 

PARAMETERS FCR NO SHOCK 
EME3=EME1 



o o o 


EME2=CHE1 
P2P1=1. 

P3P2«1. 

P3Pl=i. 

11=1 

WRITE (6,10141 II 
ILINE=ILINE*4 

IF(SHORT.EQ.l.) GO TO 4000 
WRITE(6,1039) 

DO 40 1=1,101 

WRITE (6, 10151 I,YV(I»,£M(I) ,URA T ( 1 1 , PTR AT ( I ) , PTRNS ( I) , TTRAT( 1 1 
ILINE=ILINE+1 
CALL FRIN (ILINE.Zl 
IF(ILINE-2) 38*38,40 
38 WRITEI6, 10141 II 
ILINE=ILINE+4 
WRITE(6,1039I 
40 CONTINUE 

OETERHINE UPSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 

4000 CALL FLUX (101, YY,Z,EH£1,0EL1R1 

HI1=(SHFAC-SIG1)/(TWTTE*(1.*SIG11 ) 

RETH1=R£0EL1*DDS0 
PRAN0L=.72 
PR13=PRANCL** (1./3.) 

T3ARTO=.5*TWTTE*.22*PR13+( .5-.22*PR13l/<l.*SICll 
CFLT1=.246*EXP(-1 .561*HI1) • (RET HI** (-• 268 » I / 

1 ( (TBARTO* (l.+SIGl ))**0, 79631 
WRITE (6, 1016) DSO,ODSD,BLMN,BLMON,SHFAC,HI1,CFWW1,CFLT1 
ILINE =IL1N£*6 
DS01=0S0 
00S01=00SD 

WRITE(6,1017) DEL1R,BLMI,BLM0I 
ILINE=ILINE*3 
1=1 

UTU£3A(I»=UTUE3 

FM1=BLMI 

FH01=eLM0I 

J=1 

D£L30R(J)=03120*0EL1R 

R3=BL£E0*2.*FH1 

L=1 

WRIT£(6,1036) L,BLE€01 
L=2 

WRITE(6,1036) L*BLEE02 
WRITE(6,1035) BLEE0,R3 
IF(IZ.EQ.l) RETURN 
YRS=0, 

IF(Aie.EQ.l.) GO TO 500 
IF(aL££ 01 .HQ.O. ) GO TO 500 
WSA=AES(BLE£D1) *2,*FM1 
CALL INTRP(WR,YY, HSA,YYS,161» 

CALL INTRP( VY,PHIR,YYS*PHISA,101) 

YRS=YYS*0EL1R 

WRITE (6,1037) HSA ,YRS,PHISA 
R3R=1. 
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IF(AIB.Ea.3.) R3R=1.-YRS 


AIB=1 POROUS PLATE SUCTION 
AIB=2 SLOT SUCTION 
AIB=3 SCOOP SUCTION 

5tJ0 IFCAIB.LE.1.1 BLHCM1=0. 

IFIAIB.GE.2.1 BL«0M1=-PHISA 
KB=1 

WRITE (6» 10381 KB,BLM0M1 
IFfSHORT.EQ.l.) GO TO 50 
CALL PRINaLIN£»3l 
WRITE(6,1018) 

ILINE=ILIhE+2 

ITERATION PROCEDURE FOR SOLUTION OF CONTINUITY EQUATION 
5C Rl=2. 

GAMl=(GAMMA-l.l/2% 

R2=-2. 

IF(AIB.EQ.3*) R2=R2*R3R*R3R 
P3=BLEE0»2.*FM1 
60 DEL3RP=0EL30RCJI 

IF<AIB.EQ.3*> DEL3RP=DEL3RP-YRS 
D3iP=sO£L3RP/DELlR 
REDEL3=R£DEL1»D31P 
SIG3=GAM1»EME3^EME3 
SIGMA3=SIG3/1 1.<*SIG3I 
SIGS3=SQRTCSIGNA3I 
FSIG3=(ASIN (SIGS3))/SIGS3 
VWVE3=CTWTTE»(l*+SIG3)i *♦1.76 

PX3=. 5*C1 .-UTUE3A Cl) * ( < l./BK) *ALOG f R£0£L3*ABS IUTUE3A C I )l 
1*FSIG3/VHVE3)+Cn 

CALL PRFL (UTUE3ACI) ,PX3 ,EME3 , P3P1 ♦ It YY ) 

DEL3R=DEL30RC J) 

C0N=0EL3R 

D3PR3= (DEL3R-YRS) 7R3R 
IFCAI8.EQ.3.) DEL3R=D3PR3 
IFCSHORT.EQ.O.) WRITE C6 1 1040 I R3R,03PR3 
ILlNE = ILINE<-2 

CALL FLUX(101tYYt2tEM£3tOEL3R) 

GEL3R=C0N 

FM3=BLMI 

FH03=BLM0I 

Ell J) = Rl*FMlfrR2*FM34^R3 
IFCSHORT.EQ.l.) GO TO 65 

WRITE (6,1019) Ri* R2, R3 , FMl , FM 01 , FM3, FM03 

WRITE (6,1020) I,UTUE3A( I) ,E2(I) ,J,O£L30R( J) ,E1(J) 

ILlNE=ILINE+6 

CALL FRIN(ILINE,2) 

IF(ILlNE-2) 64,64,65 

64 WRITE(6,1C18) 

ILINe=ILINE^2 

65 TEST^ABSCEIUM 
IFCTEST.LE.C.OOOCl) GO TO 80 
IF(J*G£.2) GO TO 70 

J=2 
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OEL3OR<JI=0.9*OEL3OR<l) 

GO TO 60 

70 SL 0 PE 1 = CE1C J-ll-EUJn/(0EL30RC J-l> -DEL 30RCJM 
DEL30RIJ<-1)=0EL30R<JI-E1CJ)/SL0P£1 

J=J+1 

IFCJ.EQ.SII GO TO 71 

GO ro 60 

STOP CASE NO CONVERGENCE ON J 

71 WRITEC6,1021> 

RETURN 

ITERATION PROCEDURE FOR SOLUTION OF MOMENTUM EQUATION 

dO S1=0« 

S2=0- 

S3=0. 

S4=2. 

S5=-2. 

S6=8LMOMl 

S7=0. 

T1=0# 

T2=0. 

T3=0. 

T4=0. 

IFCAIB.EQ.3.) S5=S5»R3R»R3R 

E2«I)=Sl»Tl + S2»T2^P3Pl+S3»T3<'S4»FM01frS5»FM03+S6*S7*T4 
IFCSHORT.EQ.l.) GO TO 82 
WFIT£C6tl022> S4* S5 ,T1 , T2 , T3,T4 

WRITE (6, 1020) ItUTUESAd) ♦E2CI) * J t DEL30R1J) , Elf J) 

ILINE=ILINE*'6 
CALL PRIN(ILINE»2I 
IFaLINE-2> 81»8lTa2 
dl WRITE (6*1018) 

ILINE=ILINE^2 
82 TEST =ABS(£2Cin 

IFfTEST.LE.O.OOOOl) GO TO 110 
IF(I*GE.2) GO TO 90 
1=2 

UTU£3A(I)=0*9»UTUE3 
GO TO 100 

90 SLOPE2=(£2CI-l)-E2(I))/ (UTUE3 A( I-D-UTUE3 A (I) ) 
UTUE3A(I+1)=UTU£3A(I)-E2 (II/SLOPE2 
1 = 1+1 

IF(I.EQ,51) WRITE (6, 1023) 

IF(I*£Q.51) RETURN 
100 DEL30RC1)=DEL30R(J) 

J=1 

GO TO 60 

SOLUTION TO CONTINUITY AND MC.fENTUM EQUATIONS HAS BEEN DETERMINED 
110 UTU£3=UTUE3AfI) 

DETERMINE DOWNSTREAM BOUNDARY LAYER PROFILE PROPERTIES 
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CALL PRFL (UTUE3,PX3tEHE3,P3Pl,2,YYI 
IFfSHCRT.EQ.O.) CALL PRIN C ILINE »3I 
WRITEC6tl018l 
1100 11=3 

WRITE(6»1014I II 
WRIT£(6,1C43) UTUE3»REDEL3 tPX3 
ILINE=ILINE+6 

IFCSHCRT.EQ.l.l GO TO 1200 
WRIT£I6,1039) 

DO 120 1=1,101 

WRITE (6,1 CIS) I,YYCI) ,EM(I) ,URATIII , PTRAT £ I ) , PTRNS (II , TTRATII ) 

ILINE=ILIN£+1 

CALL PRIN (ILINE, 2) 

IF(ILINE-2) 111,111,120 
111 WRITE(6, 10141 II 
ILlNE=ILINE+4 
WRITE(6, 10391 
120 CONTINUE 

DETERMINE DOWNSTREAM BOUNDARY LAYER INTEGRAL PROPERTIES 

1200 IF(AIB,£Q.3.) DEL3R=03PR3 

CALL FLUX(101,YY,2,£ME3,DEL3RI 
FUTUE3=UTUE3»FSIG3 

CFWW3=2.»FUTUE3^FUTUE3/ (T WTTE»C 1. «»SI G3 1 1 
HI3=(SHFAC-SIG3I/ (TWTTE»(1.^SIG3I ) 

RETH3=RE0EL3^00S0 

TBARTO=.5»TWTTE + * 22*PR13+ ( . 5- . 22»PR13I / (1 . ♦SIGS I 
CFLT3=*246^EXP(-1,561^HI3)^ (RET H3»* ( •. 268) I / 

1( (TBARTO^ (1.+SIG3)) ♦♦0*7963) 

WRITE (6, 1016) OSO,ODSO,BLMN,BLMON,SHFAC,HI3,CFWW3,CFLT3 

IF(AIB.EQ.3.) 0EL3R=03PR3^R3R 

WKIT£(6,iC17) DEL3R,BLMI,BLM0I 

lLIKE=ILINE+3 

D31=:DEL3R/0EL1R 

DS31=DSD^031/0SD1 

DOS 31=ODSO^D31/DOSD1 

WRITE (6, 1024) D31 , DS31, DDS31 

IFCAie.EQ.3*) WRITE(6,1041) Y RS ,D3PR3, R3R 

RETURN 

1000 FORMAT(7F10.&) 

ICll FORMAT(1HG/10X,35HSOLUTION FOR AXIALLY SYMMETRIC FLOW) 

1014 FCRMATdHO, 9X , 35H30UN0ARY LAYER PROFILE DATA STATION, 131 

1015 FORMATdH , 6X , 1 3, 6F14. 6 ) 

1016 FORMATC1HO/10X,14HDSLSTAR/OEL = , F14 . 6 ,4X , lOHMOM/DEL = , 
1F14*6/10X,28HNON-DIM£NSIONAL MASS FLUX = ,F1C.6/ 
210X,32HNON-DIMENSIONAL MOMENTUM FLUX = ,F1C.6/10X, 

315HSHAPE FACTOR = , FlC . 6, 5X , 3 OH INCOMPRE SS IBLE SHAPE FACTOR = , 

4 FIO. e/iaX,16HCF<WALL-WAK£) = , FI C .6 ,5X , 2 2HCF (LUO WE IG TILLMAN) = 
5F10*6) 

1017 FCRMATdH , 9X,8HO£L/R = ,F10 .6/1 OX, 21H3.L • MASS INTEGRAL = • 
1F10.6/13X,25H8.L. MOMENTUM INTEGRAL = ,F1D.6) 

1018 FORMAT (IHC, 9X , 42HRESULTS OF ITERATION FOR UTU£»3 AND DEL3/R) 

1019 FCRMATCIHO, 9X,5HR1 = , FIO . 6, 6X ,5 HR2 = , FI 0 . 6 , 5 X , 5HR3 = ,F10,6/ 
110X,6HFM1 = ,F10.6,5X,7HFM01 = ,F 10 . 6/ 1 OX , 6HFM3 = ,F10.6,5X, 
27HFM03 = ,Fi0.6) 

1020 FORMATdH , 9X,4HI = , I 3 , 3X , 9HUTUE»3 - , FIO . 6 , 5 X, 5H£2 = ,F10,6/ 


110X,3HJ =,I4t5X»9H0EL3/R = »F10 .6* 5X»5HE1 = »F10.61 

1021 FORHATUHO, 6X*19HN0 CONVERGENCE ON J) 

1022 FORHATCIHC, 9X,5HS4 = t FIO • 6, 5X ,5H$5 = ,F10.6/10Xt 

lEHTl = fF10*6»5X55HT2 = , FIO. 6, 5X,5HT3 = t FIO . 6t 5Xt 5HT4 = tF10.6l 

1023 FORMAKIHO, 6Xtl9HN0 CONVERGENCE ON I) 

1024 FORMAT(1HO/10X,12HO£L3/OEE1 = tFlC.6/10X, 

12aHOELSTAR3/0£LSTARl = t FlC .6/1 GX t 12HM0M3/M0M1 = tFlO.SJ 

1035 FORMAKlHGt 9X t 19H (M) BLEED/ i M)8 .L . = ,F1C*6/ 

110X*19H(M) BLEED/(M)GUCT = tFlC.6) 

1036 FORMATdHC, 9X » OH «M) BLEED t II, IIH/ t M) 8 . L. = ,F10.6, 

120H (BLEED POSITIVE INH 

1037 FORMAT(lHO/iOX,6HWSA = ,F10 .6 , 5X, 6HYRS = , FlC . 6, 5X, 8HPHISA = , 
1F10.6) 

1038 FORMATIIHO, 9X ,40 HMOMENTUM FLUX ASSOCIATED WITH B.L. BLEED, lit 
13H = ,F14.6,14H (POSITIVE INIJ 

1039 FORMATdHO, 

1 9X,1HI,8X,1HY,14X,1HM,12X,4HU/UE, 9X,6HPT/PT1, 6X, lOHPT ( NSI /PTl, 

2 6X,6HTT/TT1//I 

1040 FORHATCIHC, 10X,7HR3/R = , FIO . 6, 5X , 12H0EL3 PR/R3 = ,F10.6I 

1041 FORHATdHC, 8X,2aHSC00P CASE, SCOOP HEIGHT/R = ,F10.6,5X, 
112H0EL3PR/R3 = ,F1C . 6, 5X, 7HR3/R = ,F10.6) 

1043 FORMATdHO, 10X,9HUTUE»3 = , FIO. 6, 5X, 9H RE0EL3 = ,E14,6,5X, 

16HPX3 = ,F10.6> 

END 


SUBROUTINE SHOCK! ISH, ALP, TH ,EMX , EHY, PRAT) 

C SUBROUTINE TO CALCULATE CHANGES IN PROPERTIES ACROSS AN OBLIQUE 
C SHOCK WAVE 

DIMENSION TTRAT(2QO),PTRAT(200I,PTRNS(2CO),URATC200),EM(200) • 
DIMENSION AC(4), WR(200), PHIR(200) 


COMMON 

TTRAT 

, 

PTRAT 

, 

PTRNS 

, 

URAT 

f 


1 

EM 

, 

WR 

, 

PHIR 

, 




2 

SHORT 


BLEED 1 

, 

BLEED 2 

, 

BLEED 

, 


3 

AI8 

, 

TWTTE 

, 

GAMMA 

, 

THID 

9 


4 

BLMN 

, 

BLMON 

, 

8LMI 

, 

BLMOI 

9 


5 

OSO 


OOSD 

, 

BK 


EMEl 

, SHFAC 

f 

6 

SIGl 

, 

SIGMAl 


SIGSl 

, 

FSIGl 

, VWVEl 

9 

7 

CFHHl 


C 

* 

STCOEF 






G1=EMX»£MX 

G2=G1^G1 

IF(ISH-l) 10,10,20 
C 

C DETERMINE* SHOCK WAVE ANGLE 
C 

10 G3=CGl+2.)/Gl 
G4=(2.»G141.)/G2 

G5=(GAMMA + 1.)*(GAMMAM. I /4. + ( GA MM A-1 . ) /G1 
S0S=SINIALP)^»2. 

ACd) = (S0S-l.)/G2 
AC(2I=G4+G5*S0S 
AC(3)=-G3-GAMMA»SDS 
AC(4)=1. 

C 

C FIND ROOTS OF CUBIC EQUATION IN <S IN ( THET A ) I **2 
C 
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CALL CUBIC(AC,Z) 

ROOr=SQRT(ZI 
IF(R00T-1,I 12»12,11 
C 

C ERROR EXIT IF SINCTHETAI GREATER THAI# r.ff 

C 

11 WRIT£I6,1C10I ROOT 
TH=ASINCROOn 
RETURN 

12 TH=ASIN(ROOT> 

C 

C DETERMINE PRESSURE RATIO ACROSS SHOCK ANO DOWNSTREAM MACH NO. 
C 

20 GAM1=GAMMA^*1. 

GAM2=GAMMA-1. 

GAM3=GAM1*GAM1 

2=SIN(TH>*»2. 

PRAT= C2.»GAMMA*G1*Z-GAM2)/GAM1 
ANUM=GAM3»G2»Z-4.*CG1*Z-1.)*CGAMMA^G1^Z+1.I 
DENO= f2.»GAMMA»Gl*Z-GAM2)»CGAM2*Gl*Z«-2.l 
EMY = SORT< ANUM/DENOI 
RETURN 

1010 FCRMAT(1HC,10X,10HSIN«TH) = *F10.6I 
END 


SUBROUTINE FLUX tK . Y t INDIC t EME »OEL R) 

C SUBROUTINE TO CALCULATE MASS AND MOMENTUM FLUX OF B.L. 

C ALSO CALCULATES DISPLACEMENT AND MOMENTUM THICKNESSES 
C INDIC=1 TWO-DIMENSIONAL FLOW 

C =2 AXIALLY SYMMETRIC FLOW ON COWL 

C =3 AXIALLY SYMMETRIC FLOW ON CENTERBOOY 

DIMENSION TTRATfZOO) .PTRATfZQO) »PTRNS(230) , URATI2C 0 ) ,EM C 200) 
DIMENSION YC2aO),YY{20C) ,BLMR<200) ,BLMOR<2PO) tBLMIR<200) 
DIMENSION PMOIR(2GO),YR(2flO),WR(200)* PHIRC200) 


COMMON 

TTRAT 

* 

PTRAT 

. PTRNS 

t 

URAT 

• 


1 

EM 

• 

WR 

• PHIR 

« 




2 

SHORT 

« 

BLEED 1 

« BLEED 2 

V 

BLEED 

f 


3 

A IB 


TWTTE 

« GAMMA 

f 

THIO 

« 


4 

BLMN 


8LM0N 

« BLMI 

» 

BLMOI 

t 


5 

OSD 

t 

OOSD 

• BK 


EMEl 

« 

SHFAC 

6 

SIGl 

• 

SIGMAl 

t SIGSl 

• 

FSIGl 

t 

VWVEl 

7 

CFWHl 


C 

. STCOEF 





DO 100' 
PRAT=1. 

I ItK 









TTOT=l.+CGAMMA-l. )*EMCI )*EMII)/2. 
TTOTE=l.+(GAMMA-l.)^£ME»EME/2. 
TOTE=TTOTE»TTRAT« D/TTOT 
RHRAT=PRAT/TOT£ 

GO TO (13*2Q«30)«INDIC 
C 

C TWO-CIMENSIONAL FLOW 
C 

10 BLMR(I)=RHRAT»URATII) 

BLHOR(I)=0LMR(I)*URAT Cl) 

IFCURATCI) .LE.O.) BLMOR ( I ) --BLMOR (I) 



ooo ooo noo ooo 


YY(II=Y«I) 

GO TO lOD 

AXIALLY SYMMETRIC FL CH ON COWL 

20 BLMR(II=RHRAT^URATCI) 

BLMORdlsBLMRdl^URATai 

IFIURAT(I) .LE.O.) 8LH0R (I)=-BLM0Rf II 

YYdl=Y<I) 

YRm=YYII)*0ELR 
BLHlRCI)=BLMRd)»fl.-YRd)l 
BM0IR(I) = BLM0R«I) *(l^-VR(in 
GO TO 100 

AXIALLY SYMMETRIC FLOW ON CENTERBODY 

30 GO TO 100 
100 CONTINUE 

GO TO <110tl20,130ldNDIC 

THO-CIHENSIONAL FLOW 

110 DO 115 1=1, K 

call INTEG<I,YY,BLMR,AREA1I 
CALL INTEG(I,YY,BLM0R,AREA2I 
WRd)=AREAl 
PHIR(I)=AR£A2 
115 CONTINUE 
BLHN=AREA1 
BLM0N=AREA2 
DSD=1.-BLMN 
DOSO=ELMN-BLMON 
SHFAC=DSD/OOSO 
RETURN 

AXIALLY SYMMETRIC FLOW ON COWL 

120 00 125 1=1, K 

CALL INTEG(I,YR,8LMIR,AREASI 
CALL INTEGCI,YR,BM0IR,AREA61 
WR(Ii=AREA5*2. 

PHIRdl=AREA6*2. 

125 CONTINUE 
BLMI=AREA5 
BLM0I=AR£A6 

BLMN=2**8LMI/(2.»rELR -OELR^OELRI 
BLM0N=2*^ELM0I/<2.*0ELR -OELR^OELRI 

DSD=(l.-SCRr(i.-2,»0£L Pf-OELR^DELR *2 .-^BLMII ) / DELR 
0aS0=(l.-SQRT(l,-2.^3LMI«*2,»BLH0I 1 1 /DEL R 
SHFAC=OSO/ODSO 
130 RETURN 
ENO 


SUBROUTINE PRFL (UTUEST , PX , EME,PKP 1, lOPT, YY ) 

C SUBROUTINE TO CALCULATE DISTRIBUTIONS OF PROPERTIES 
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C FOR BOUNDARY LAYER WITH WALL-WAKE tfELOCITY PROFILE 

DIMENSION TTRATC2C0) •PTRAT<200) ,PTRNS(2»50I # URATC200 l»EM (20fl) 
DIMENSION YY(200I^WRC200)» PHIRI200) 


COMMON 

TTRAT 


PTRAT 

9 

PTRNS 

9 

URAT 

9 



1 

EM 


WR 

9 

PHIR 

9 





2 

SHORT 

• 

8LEE01 

9 

BLEE02 

9 

BLEED 




3 

AIB 

• 

THTTE 

9 

GAMMA 

9 

THIO 

9 



4 

3LMN 


BLHON 

9 

BLMI 

9 

BLMOI 

9 



5 

DSD 

f 

OOSO 

9 

BK 

9 

EMEl 

9 

SHFAC 

• 

6 

SIGl 

f 

SIGMAl 

9 

SIGSl 

9 

FSIGl 

9 

VHVEl 

9 

7 

CFWHl 


C 

9 
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COMMON /UWU/AIOO 
PI=3. 1415927 
EXP2=GAMMA/CGAMMA-1.I 
EXP3=1. /(GAMMA-1. I 
GAMl=(GAMMA-l.)/2. 

GAMZ^GAMMA^-l. 

GAM3=GAMMA-1. 

G1=EM£*£ME 

SIGMA=GAM1»G1/(1.+GAM1*G1I 

SIGG1=SQRT(5IGMAI 

SIGG2=1./SIGG1 

SICG3=ASIN(SIGG1I 

URAT(1)=0. 

TTRATCH = TWTTE 
EM(1) = 0. 

YY(1)=0. 

PTRAT ( 1) = Cl •/( 1 •♦GAM1»EM£1*EME1 IJ ♦♦EXP2»PKP1 
PTRNSC1)=PTRATC1) 

DO 100 1=2*101 
A 1=1-1 

YYf II =AI/100. 

YAA=SCRT(l.-YY(II^»AIDOI 

ALG=ALOG( YY(I) H'2 .^ ( Y AA-ALOG C 1. +Y A All /A lOD 
URAT(II=SIGG2^SINCSIGG3-SIGG3»PX* ( 1 . *0 OS (PI »YY ( 11)1 ♦ 
1(1./BKI»UTUEST*SIGG3»ALG) 

TTRAT(I)=TWTTE+<l.-TWTTE)*A8S(URATfI)| 

U2=URAT(I)^URAT(I) 

£M(I) =SQRT (U2/ ( (1 ./Gl+G AM 1 I ♦TTR AT ( II -G AM1*U2) I 
IF(URAT(II.LE.O.I EM ( II =-l. ♦EMC II 
IFCIOPT-1) 100*100*20 
C 

C CALCULATION OF TOTAL PRESSURE DOWNSTREAM OF NORMAL SHOCK 
C 

20 PTRAT C 11 = C Cl. ♦GAMl^EM (I) ♦EM <11 1 /C 1 . ♦GAM1*EME1^EM£1 1 1 ♦♦ EXP2»PKP1 
IFCEMIII-I.I 40*40*50 
40 PTRNSCIJ=PTRATCII 
GO TO 100 

50 PTRNS CI)= (GAM2»EM (I) ♦EMCII/2./C 1. +GAM1^ EME1»EME 11) ♦♦EXPE^ C GAM2/ 
1(2.»GAMMA»EMCI)^EMCII-GAM3II»^£XP3^PKP1 
ICO CONTINUE 
RETURN 
END 


SUBROUTINE INTRPC X* Y* XX* YY* Nl 
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C LINEAR IMERPOLATION 

DIMENSION XC200NY(200) 

DC 30 1=1 tN 
IF(XX-X(I)J 10,20»30 
10 K=I-1 

YY=YCK)*(Y (I )-Y(K)J »4XX-XIK11/CX<I J-XCKI) 
RETURN 
20 YY=Y4I) 

RETURN 
30 CONTINUE 
RETURN 
ENO 


SUBROUTINE CUBICCC^ZI 

SUBROUTINE TO DETERMINE ROOTS OF A CUBIC EQUATION 
NASA» AMES RESEARCH CENTER, MOFFETT FIELD, CALIF. 
DIMENSION C(<*),VLSTf 7) 

ACOS(X) = ATAN<SQRT Cl . 0-X*»2l/XI 
P=-C(3)»*2/3.0 ♦ CC2) 

Q=2.0^C(31^»3/27. n - C ( 2) »C C3 )/3. 0 ♦ CCl) 

RSQ = -C.5»Q/ SQRT(-F*^3/27.0I 

IF (ABS(RSQ) .GT. 1.0) RSQ=S IGN ( 1. 0, RSQI 
PHI=ACOSI RSQ) 

TEM=2 .0^SGRT(-P/3.9» 

Xl= TEM»COS(PHI/3.0) 

X2= TEM»COS(PHI/3.0 ♦ 2.99439510) 

X3= T£M»COS(PHI/5.0 ♦ 4,ld879C20) 

IF (X2-X3) 153,150,160 
150 Yi=AMAXl(Xl,X2i 
Y1=AHIN1C Y1,X3) 

GO TO 175 

160 Yl=AMINlf X1,X2) 

Y1=AMAX1( Y1,X3I 
175 Yl=Yl-CC3)/3.0 
8 Z=Y1 
RETURN 
ENO 


SUBROUTINE INTEG ( K, Y , Z, ARE Al 
INTEGRATION USING SIMPSON^S RULE 
DIMENSION Y(230), Z(20C) 
IF(K.GE,5) GO TO 1 
IFCK.EQ.l) GO TO 21 
IFCfC.EQ.Z) GO TO 22 
IFIK.EQ.3) GO TO 23 
IFCK.EQ.M) GO TO 24 
1 AK=K 
BK=AK/2. 

KK=BK 

CK=KK 

IF(BK-CK) 4,2,4 
K IS EVEN 
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c 

2 N=K-1 
GO TO 5 

K IS ODD 

4 N=K 

5 O0Q=0. 

EVEN = 0. 

J=N-3 

DO 10 1=2, J ,2 
EVEN = EVEN ♦ Z fll 
000 = 000 ♦ Z(I+11 
10 CONTINUE 

AREA = fY(2)-Y(l>)/3**(Z(l)+ZCNI+4.»IEVEN*Z(N-l)l+2*»000l 
IF(BK-CK> 14,12,14 

K IS EVEN 

12 AREA = AREA ^ I Y( Kl-Y (K-1) I ♦ < Z« K> ♦ZC K-1 >) /2* 

RETURN 

K IS ODD 

14 RETURN 

21 AREA=0* 

RETURN 

22 AREA=<Y(2)-Y(1) )^(Z(2)+Z<1) 1/2* 

RETURN 

23 AREA= (Y{2)-Y(l) )^(Z(3) +4 • »Z ( 21 ♦Z C II ) /3* 

RETURN 

24 AREA=(Y(2)-Y(1I !♦ (CZC4) l3)l/2,+ CZ ( 3 ) *4. ♦Z (2 ) + Z ( 11 l/3#l 
RETURN 
ENO 


SUBROUTINE PRIN ( ILINE , I NO I 

SUBROUTINE TO COUNT LINES OF OUTPUT ON PAGE, CHANGE PAGE WHEN 
NECESSARY AND WRITE TITLE ON EACH PAGE 
DIMENSION TITLE(121,TITL1(12I 
GO TO (10,2J,30I,INO 

FIRST PAGE TITLE 

10 READ(5, 10001 CTITLEIJI ,J=1, 121 
REAO(5,10 00J (TITLKJ) ,J=1,12) 

ILINE=2 

WRITE(6,101P) (TITLECJ), J=l,12l 
WRITEC6tl012) (TITLK J) * J=l,12) 

IPAGE=1 

WRITE(6. 10111 IPAGE 
RETURN 

TEST FOR END OF PAGE 
20 IFdLINE-48121,22,22 
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21 RETURN 
C 

C CHANGES PAGE 

C 

22 ILINE=2 
IPAGE=IPAGE+1 

C 

C PAGE TITLE 

C 

WRITE«6,lfllO) (TITLE(J>, J=l»12) 
WRITE <&♦ 1C 12) (TITLiCJ) f J=l,12) 
WRIT£(6»1011) IPAGE 
RETURN 
C 

C PROGRAMMED PAGE CHANGE 

C 

30 IPAGE=IPAGE+1 

HRITE(6,1Q1Q) < TIT LE ( J) » J=l, 12) 

WRITE(6,lul2) (TITLKJ)» J=l*12l 

WRITE(6»1011) IPAGE 

ILlNE-2 

RETURN 

1000 F0RMAT(12A6) 

1010 FORMAT(1H1,30X,12A6) 

1011 FCRMATdH ,100X»5HPAGE »I3I 

1012 FCRMATCIH t30X,12A6) 

ENO 


TABLE 3-B: INPUT TO PROGRAM ANAL 


1C OEGPEt CONE 
CHEN CHH SUN 

NC ALPHA=1. 


2,82 0 

.0446 

71200. 

1. 

O.C 0 

.0 

0. C 

C.4 

22.7 

0.15 2 

.0 

0.0 


50. 

22.751213 

2.902980 

2.684648 

1.229258 

22.711516 

3.182413 

2.673301 

1.250928 

22 .091668 

3.3C6264 

2.668448 

1.260 316 

21.671819 

3.423d2j 

2.663971 

1.269043 

21.651971 

3.534813 

2.659794 

1.277 243 


0.0 
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21,47ei22 

21.212274 

20.992425 

20.772577 

20.552728 

20.113031 

19.893183 

19.673334 

19.453486 

19.233637 

18.793940 

18.574092 

18.354243 

18.134395 

17.914546 

17.474849 

17.255C01 

17.035152 

16.815304 

16.595455 

16.155758 

15.935910 

15.716061 

15.496213 

15.276364 

14.836667 

14.616819 

14.396970 

14.177122 

13.957273 

13.517576 

13.297728 

13.077879 

12.858031 

12.638182 

12.198465 

11.978637 

11.758788 

11.538940 

11.319391 

10.879394 

10.659546 

10.439697 

10.219849 

lO.OCClOO 


3.642199 
3.746224 
3.847539 
3.946660 
4.344C05 
4.234706 
4.328607 
4.421846 
4.514618 
4.607097 
4.791800 
4.884306 
4.977087 
5.07C264 
5.163953 
5.353313 
5.449201 
5.546336 
5.643924 
5.742972 
5.944979 
6. 048160 
6.152942 
6.259446 
6.367793 
6.590529 
6.705188 
6.322235 
6.941821 
7.064109 
7.317486 
7.448952 
7.583874 
7.722475 
7.864991 
8.162819 
8.318704 
8.479660 
8.64604J 
8.818227 
9.181742 
9.374038 
9.574088 
9.782512 
10.3C0Q0G 


2.655863 

2.652139 

2.648592 

2.645198 

2.641939 

2.635768 

2.632832 

2.629984 

2.627216 

2.624522 

2.619333 

2.616830 

2.614383 

2.611988 

2.609643 

2.605093 

2.6G2885 

2.6G0718 

2.598592 

2.596506 

2.592451 

2.590480 

2.538546 

2.586650 

2.584792 

2.581190 

2.579447 

2.577744 

2.576082 

2.574463 

2.571360 

2.569879 

2.56845C 

2.567074 

2.565756 

2.5633C6 

2.562184 

2.561137 

2.560173 

2.559297 

2.557845 

2.557286 

2.556860 

2.556574 

2.556446 


1.285010 
1.292413 
1.299505 
1.306328 
1.312916 
1.325487 
1.331511 
1.337 382 
1.343113 
1.348717 
1.359576 
1.364847 
1.370022 
1.375105 
1.380101 
1.3B9849 
1.394607 
1.399290 
1.4039D1 
1.4C8441 
1.417312 
1.421644 
1.425907 
1.430100 
1.434222 
1.442248 
1.446140 
1.449969 
1.453707 
1.457 359 
1.464386 
1.467 751 
1.471007 
1.474147 
1.477164 
1.482786 
1.485369 
1.487782 
1.490 010 
1.492 035 
1.495400 
1.496691 
1.497 685 
1.498350 
1.498649 



TABLE 4 


INPUT TO METHOD OF CHARACTERISTICS PROGRAM 


rtACH 

2.62 10 

□£6RE£ CONE 


2.05 




0 

47 1 

0 0 0 


2.62 

1.4 

23.9999 22.78 

2.5 

0.QC0C31 

U .001 

0.00001 .01 


3.5 

0.001 

0.001 0.0 

5.0 

2. 

2. 



0.0 

5.0 



0.0 




10. 

10. 



l.C 

0.99 

5.0 



0.0 

Q.O 



0.0 

1. 

2.82 


2 

3 2 

1 2 22 



1.6350000 

•9900000 

Q.oooao 2 

• 82 

1.0 

1.635C00G 

.9800000 

0.009DC0 2 

• 82 

1.00 

1.635Q00U 

.9700000 

0.9000 2 

• 82 

l.QO 

1.635Q000 

.96CC0Q0 

O.rOQOO ^ 

• 82 

1. DO 

1.8359000 

.95C0030 

0.9CC093 2 

• 82 

1.00 

1. 6353000 

.9490900 

O.^'OCOQ 2 

• 82 

l.CQ 

1. 6350000 

.9300000 

0.C030C032.62 

1. 

1.8350000 

.9209030 

O.OOOCO 2 

.82 

l.CO 

1.835C000 

.91COOOQ 

O.OOGOO 2 

.82 

l.OQ 

1. 8350030 

.9000000 

O.GGCOO 2 

.62 

1.00 

1.83500C0 

.8900000 

0.00000 2 

• 82 

1.00 

1.8359000 

•8900002 

O.COOOO 2 

• 82 

1.00 

1. 835C000 

.8790000 

O.nOQOC 2 

.82 

1.00 

1.8350300 

•8600003 

O.OOOQCO 2 

.82 

l.CO 

1.8350000 

•8500000 

0.0000 1 

.82 

l.OQ 

1.835C0QC 

.84CQOOO 

O.QQOQO 2 

• 82 

1. 00 

1.835COOO 

.8300000 

0.0090C 2 

.82 

1.00 

1.8350^00 

.8209000 

0.0000 2 

.82 

1.00 

1.835G903 

•dlCCCDO 

O.C09000 2 

.82 

1.00 

1.835C000 

.8300003 

9.0QC30QO 

2.82C0CCn 

1.3030000 

1.8350000 

.79COOOG 

0.PQ30GQC 

2.82CQ uCC 

l.OOOGOOQ 

1.8350900 

.7890000 

0. .0Q300C0 

2.8200000 

l.OQQOOOO 

1 3 

3 1 

1 27 



1.635000 

♦ 323560 

10 .000 000 

2.556446 

.999057 

1.835000 

♦345365 

9. 374^38 

2.557266 

.999057 

1.835000 

.367305 

8.818227 

2.5 59 29 7 

.999057 

1.335Q00 

.389326 

8.318704 

2.562184 

.999057 

1.835C00 

.411455 

7.864991 

2.565 756 

.999057 

1.835C00 

.433699 

7.448952 

2.569879 

♦999057 

1.835000 

.456364 

7.0641C9 

2.574463 

♦999057 

1.835C00 

.478557 

6.705188 

2.579447 

.999057 

1.835C00 

•501186 

6.367793 

2.584792 

♦999057 

1.835309 

.523957 

6 .Q48160 

2.590480 

. 999057 

1.835C00 

.546879 

5.742972 

2.596506 

♦999057 
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1.83500Q 

1*835000 

1.835000 

1.835000 

1.835G00 

1.835003 

1.835000 

1.835000 

1.835000 

1.835300 

1.835000 

1.835000 

1.835000 

1.835000 

1.835000 

1.835000 


.569958 
.593203 
•616622 
.628399 
.640223 
.652094 
•664014 
.675983 
•688004 
.700C7d 
. 712202 
.724382 
.736617 
.748909 
.761259 
.773668 


5.449201 

5.163953 

4.8843Q6 

4.745611 

4.607097 

4.468279 

4.328607 

4.187440 

4.044CG5 

3.897346 

3.746224 

3.588973 

3.423220 

3.245340 

3.049241 

2.823021 


2.602885 
2.6C9643 
2.616630 
2.620607 
2.624 522 
2.628590 
2.632832 
2.637271 
2.641939 
2.646877 
2.652139 
2.657800 
2.663971 
2.670822 
2.678 641 
2.687994 


.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 

.999057 
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TABLE 5-A 


PROGRAM BLGRN 


PROGRAM BLGRN (INPUT, TAPE 5=INPUT,0UT»»aT*T«l»f6=0UTPUTl 
INPUT FORMAT 10F7.0 EXCEPT CARD 1 
CARO(S) 

1 title, columns 1-72, HOLLERITH 

2 COL.i-7 AJ AJ-C. FOR 2-D FLOW 

A4«l. FOR AXISYMMETRIC FLOW 
6-lA TO TOTAL TEMPERATURE (R) 

15-21 TW WALL TEMPERATURE (R) 

22-26 REYNOF FREE STREAM REYNOLDS NO. PER FT /lUOOOO. 

29-35 SIGMA PRANDTL NUMBER (-.72) 

36-A2 CF COEFFICIENT OF SKIN FRICTION 

A3-A9 DEL BOUNDARY LAYER THICKNESS (IN.) 

3 COL. 1-7 XO VALUE OF X FOR STARTING POINT (HAY EQUAL TO ZERO) 

e-lA DXMN MINIMUM PERMISSIBLE STEP SIZE (IN.) 

15-21 DXMX MAXIMUM PERMISSIBLE STEP SIZE (IN.) 

22-28 ANP FREQUENCY OF PRINTOUT AS MULTIPLE OE STEP SIZE 

29-35 ANXP NUMBER OF POINTS AT WHICH EXTRA PRINT IS REQUIRED 

4 COL. 1-7 XEND VALUE OF X FOR ENDUN6 POINT 

8-14 ASEA PARAMETER IN ( l.-( Y /DEL ) ♦♦ASEA ) (-1. IS RECOMHENDE 

5 COL. 1-7 ANME NUMBER OF VALUES GIVEN IN MACH NUMBER ARRAY(BETWEE 

4. AND 200.) 

8-14 ANSME NUMBER OF SMOOTHING PASSES ON MACH NO. DATA 

6 COL. 1-70 AHE VALLES OF MACH NUMBER ARE INPUT, MORE THAN ONE CAR 

CAN BE USED IF NECESSARY 

7 COL. 1-70 X CENTERLINE DISTANCE ( IN. ) AT WHICH MACH NUMBERS 

)RE INPUT 

e COL. 1-7 AMR NUMBER OF VALUES GIVEN IN RADIUS ARRAY (AT LEAST 4 

6-14 ANSR NUMBER OF SMOOTHING PASSES ON INPUT DATA 

9 COL. 1-70 R VALUES OF R ARE INPUT. MORE THAN ONE CARD CAN BE 

USED IF NECESSARY 

10 COL. 1-70 X CENTERLINE DISTANCE (IN) AT WHICH R VALUE 

ARE INPUT 

C0MM0N/XCH6E/SIGN1,P,FF,T0RR 

C0MMCN/H0NE/H1,H2 

C0MMCh/ANS/THTA,DELSTA,GASCT,BMU,GAM21 

COMMON/ APRR/ KR, NR, XIR ( 200), AIR (200), DIR ( 200), XCR( 200) 

C0MM0N/BBDD/AJ,AAM£(200),AAR (200) 

C0HM0N/C0NST/R,S,6,PI,GW,T0W,B,E91,E92,H,AME,DME,PT3,REY3,TE,IRY 
COMMON/ INTEG/CF, DEL, DCF, DDL, OX, DXMN, DXMX, X,XENO, MEND 
CtiMMON/RFIN/TO,TW,SIGMA,GAMMA,REYNO, ASEA 

COMMON/ARRME/KME,NME,XIME(2CO),AIME(200),DIME(200),XCME(200) 

925 F0RMAT(10F7.0) 

904 F0RMAT(1X,13F10.6) 

903 F0RMAT(//7X,1HX,7X,5HTHETA,5X, 8HDISP.TH.,3H H,9X,2HH2, 8X, IHP, 9X, 1 
1HF,7X,9HR*TH/1000,3X,2HME,6X,5HDELTA,5X,2HCF,6X,7HDDEL/0X,4X,6HDCF 
2/DX//) 

902 F0RMAT(//17H INITIAL DELTA -,F9.6,7H INCHES, 5X, 13H INITIAL CF «,F 
19.6) 

901 FORMAK* TO • *,F7.2/ 
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1 

♦ 

TW - 

♦#F7,2/ 

z 


SIGMA « 

♦,F6.3/ 

3 

♦ 

GAMMA * 

♦^F6,3/ 

4 

♦ 

PSIA - 

*^F7,3/ 

4 

♦ 

REYNC » 

♦^F9.2// 

5 

♦ 

XO * 

♦,F7.2/ 

6 

♦ 

DXMIN * 

♦/F7.4/ 

7 

* 

DXMAX « 

♦#F7,4/ 

8 

* 

XBND ■ 

+,F7.2/ 

9 

* 

NPRINT ■ 

♦^13) 


40 F0RMAT(8A10) 

41 F0RMAT( //#6A10//) 

49 FORMATdHl) 

12344 F0RMAT(16H PROGRAM RAN AT I2jlH.d2«3H. «2A10///) 

200 FORMAT(1H115X»100HTEM235 • A COMPUTER PROGRAM TO PREDICT THE OEVEL 
lOPMENT OF A COMPRESSIBLE TURBULENT BOUNDARY LAYER ) 

DIMENSION TITLE(8) 

CALL CL0CK(IHR»MIN/ ISEC) 

CALL DATE(A7#B7> 

WRITE(6»200) 

WRITE(6^12344) IHR»M1N>A7>B7 

5 CONTINUE 

C READ HEADING 

READ(5»40) (TITLE(I)»I«1* 8) 

IF(E0F>5)6,7 

6 CALL EXIT 

7 CONTINUE 

READ (5» 925 ) A J«TO> TW« REYNOF* SIGMA» CF*CEL 
IRY«0 

T0LRER«.00C08 
TORR«. 0000005 
GAMMA»1.4 
REYN0»REYN0F/12. 

REYNO-REYNO+100000. 

READ(5>92S)X0»DXMN«DXMX»ANP« ANXP 
NPRINT«ANP+0.1 

C NXP IS NUMBER OF POINT AT WHICH EXTRA OUTPUT IS REQUIRED. 

NXP-ANXP+0,1 
READ(5.925) XEND^ASEA 
NME • 0 

WRITEC6#41) (TITLE(I>#I«1»8) 

C READ AND PRINT MACH NUMBER DISTRIBUTION 
KME«NME+1 

READ(5*925) ANME.ANSME 
NME«KME+IFIX(ANME-0.9> 

IF(ANSME)230.225.225 

225 NSME-ANSME+0.1 
READ{5.925)( AIME(I).I«KME.NHE) 

READ (5.925) (XIME(I)»I*KME.NME> 

DO 3368 I«KHE.NME 

3388 AAME(I)-AIME(I) 

C CHECK THAT X IS INCREASING 
KP-KME+1 
DO 228 I«KP.NME 

IF (XIME(I)-XlME(I-l) ) 226.226.228 

226 WRITE (6.934) XIME ( I-l ). XIHE ( I ) 
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IE— 1 

228 CONTINUE 

CALL SHOOTH(NSHE>KHE>NHE*XIME»AIHE>DIHE) 

NR«0 

C READ AND PRINT CONTOUR 
AOO KR»NR+1 

READ (5>92S) ANR>ANSR 
NR=KR+IFIX(ANR-0.9) 

IF(ANSR) 410>420»A20 

420 NSR*ANSR+0.1 
READ(5>925) ( AIR(I)#I«KR>NR) 

DO 3399 I*KR»NR 

3399 AAR(I)«AIRm 

READ (5»925) tXIR(I)>I«KR»NR) 

C K THAT X IS INCREASING 
KP«KR+1 

DO 422 I*KP»NR 

IF (XIR(I)-XIR(I-D) 421>421>422 

421 WRITE <6»934) XIR ( I-l ) » XIR ( I ) 

IE— 1 

422 CONTINUE 

C SMOOTH RADII AND PRINT SMOOTHED VALUES 
CALL SMOOTH (NSR>KR>NR>XIR>A1R»0IR) 

IF(KR.NE.l) 60 TO 425 
XCR(1)*XIR(I) 

KR«KR+1 

425 DO 430 I»KR»NR 

XCR(I)«XCR(I-1)+SQRT( (XIR(I)-XIR(I-1))**2+(AIR(I )-AIR(I-l) )**2) 
430 CONTINUE 

C OBTAIN ME AS A FUNCTION OF ARC LENGTH 
DO 440 I»KME>NME 

CALL PI(KR,NR/XIR,XCR»XIME(I)»XCME(in 
440 CONTINUE 

CALL PI(KR»NR»XIR»XCR»XO,X) 

C CONSTANT 

C SQRT(2.)=1. 414214 

R-53.35 
G-32.174 
S»0.5fGAMMA-0.5 
Pl«3. 141592 
6W«TW/T0 
TOW-TO/TW 
B-TOW-1. 

E91 «1./0.915 
E92«E91-1. 

CV«4290. 

T015»2.27*T0**1.5 

BMU=TO15*l.E-6/(T0+198.6) 

6ASCi3» 6AMMA/(2.*S*CV*T0) 

GASQT«SQRT(6AS03) 

GAM21*(GAMMA-2.)/(6AMMA-l. ) 

DX«DXMN 

MEND«0 

NPC«1 

CALL DERIV 

PT3=REYN0*12./(REY3*AME) 

PTT-PT3/144. 
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RRRT»REYND*THT A/1000. 

WRITE (6^901) T0#TW#SIGMA^6AMMA^PTT^REYN0,X0#DXMN#DXMX#XEND/NPRTNT 

WR1TE(6^902)DEL^CF 

WR1TE(6,935) ANME#ANSHE 

935 FORMAT (+ONME-+# F7 .2# * NSME « *#F7.2) 

WRITE(6#936) 

936 FORMAT C* ME VALUES*) 

WRITE(6^937) ( AAME ( I) # I-KME^ NME ) 

937 FORMAT (1X,10F10.3) 

WRITE(6,938) 

938 FORMATC* CORRESPONDING X VALUES*) 

WRITE(6>937) (XIME(I), I-KME^NME) 

93A FORMAT(*OTWO SUCCESSIVE X FOUND WHICH ARE NOT INCREASING *y 
1 2F12.A) 

C SMOOTH MACH NUMBER AND PRINT SMOOTH VALUEW 
230 WRITE (6,939) 

939 FORMAT (♦ SMOOTHED MACH NUMBERS*/ 

1 * X ME DME*) 

WRITE (6,9^0)( (XIME(I),AIME(a),DIMECI))#I*KM£#NHE) 

9A0 FORMAT (1X,3F10.3) 

C 

WRITE(6,955) ANR,ANSR 

955 FORMAT (*0NR **,F7,2,* NSR - *,F7.2) 

WRITE(6,956) 

956 FORMATC* R VALUES*) 

WRITE (6,937) ( AAR ( I ) , I«KR, NR) 

WRITE(6,938) 

WRITE (6,937)(XIR(I),I»KR,NR) 

^10 WRITE (6,959) 

959 FORMAT (* SMOOTHED RADII*/ 

1 ♦ X R DR*) 

WRITE (6,9^0)( (XIR(I), AIR(I),DIR(I )),I-KR,NR) 

WRITE(6,903) 

C 

WRITE(6,90A)X,THTA,DELSTA,H,H2,P,FF,RPRT,AME,DEL,CF,DDL,0CF 

IRY»2 

^2 CONTINUE 

10 CALL MERSON 

11 IF (MENO.EQ.l) 60 TO 16 

IF (NPC.EQ.NPRINT) GO TO 16 
NPC*NPC+1 
GO TC 10 
C RESET COUNTER 
16 NPC*1 
C 

RRRT-REYN0*THTA/100O. 

WRITE(6,90A)X,THTA,DELSTA,H,H2,P,FF,RRRT,AME,DEL,CF,DDL,DCF 
IF( (•CF.GT.0.).AND. (CF.LT.TOLRER)) GO TO A3 
IF((CF.LT.0.).AND.(DX.LT..CC0065)) GO TO A3 
IF (MENO.EQ.O) GO TO 10 
GO TO 88 
A3 CF*-CF 

DXABS*2.*ABSCCF)/ABS(DCF) 

X-X+DXABS 

DEL*DEL+DDL*DXABS 

TOLRER=ABS(CF) 

CALL DERIV 



RRRT«P€YN0*THTA/100O. 

WR1TE(6»904)X» THTA,DELSTA»H#H2»P»FF» RRRT»AME»DEL>CF, DDL»OCF 
60 TO 42 
86 CONTINUE 
NRITE(6>49) 

GO TO 5 
END 


SUBROUTINE NERSON 
C MERSON INTEGRATION 

C V IS THE VARIABLE ARRAY (BOTH INPUT AND OUTPUT)# D IS THE DERIVATIVE 
C ARRAY. 

C OX IS CURRENT STEP LENGTH# DXHN AND DXNX NININUN AND HAXIMUH VALUES OF 
C DX. X IS INDEPENDENT VARIABLE# AND XENO ITS UPPER LIMIT. 

C 

C0MM0N/C0NST/R#S#G#PI#GW#T0W#B#E91#E92,H» AME#DHE#PT3#REY3#TE#IRY 
COMMON /INTE6/V(2)#0{2)#DX#0XMN#DXMX#X#XEN0#MEN0 
DIMENSION U(2)#VE(2)#DA(2)#DC(2)#DD(2)#DE(2) 

C BU#BL ARE UPPER AND LOWER BOUNDS ON TRUNCATION ERROR. 

DATA BU#BL/0.0005#0. 00005/ 

C 

C MARK IS SET NON-ZEP.O TO PREVENT INCREASING THE STEP LENGTH IF EITHER 
C XENO IS REACHED# OR THE STEP LENGTH HAS JUST BEEN DECREASED. 

MARK*0 

IF (X40X-XEND) 5#7#7 
C PREVENT OVERSHOOTING X END. 

7 HARK*1 
DX-XENO-X 
MEND-1 

C SAVE POINT A. 

5 DO 10 I*l#2 
10 U(I)-V(I) 

C 

15 DX2-0X/2. 

DX3-DX/3. 

DX6-DX/6. 

DX8-0X/8. 

C FIND POINT B. 

CALL DERIV 
DO 20 I-l#2 
DA(I)«D(I) 

20 V(I)«U(I)+D(I)*DX3 
C FIND POINT C. 

X-X+DX3 
CALL DERIV 

V(1)«(D(1)+DA(1))*DX6+U(1) 

V(2)=(D(2)+0A(2) )*DX64U(2) 

C FIND POINT 0. 

CALL DERIV 
X-X+DX6 
DO 40 I«l#2 
DC(I)«D(I) 

40 V(I)-U(I)+0X8*(0A(I)+3.*D(I)) 

C FIND POINT E. 

CALL DERIV 
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X-X+DX2 
DO 50 I»1^2 
DD(l)«D(n 

VCI)«U(I)+DX2*(DA(I)-3,*DC( J) + ^.*D(D) 

50 VE(I)»V(I) 

C FIND POINT F. 

CALI DERIV 
MINC«0 
HDEC-0 
DO 60 1*1^2 

V(I)=U(I)+DX6+(D4(I) + 4.*DD(I) + D(I) ) 

IF (VEtl)) 51#60ii51 

C TRUNC IS A MEASURE OF THE TRUNCATION ERROR. 

51 TRUNC-ABS{l.-V(I)/VE(in 
IF (TRUNC-BU) 54.5^,52 

52 MDEC=MDEC+1 

54 IF (TRUNC-BL) 56/56/60 
56 MINC-MINC+1 
60 CONTINUE 

IF (MOEC.GT.OI GO TO 70 
IF (MINC.EQ.2.AND.MARK.EQ.0) GO TO 80 
65 RETURN 

C IF EITHER TRUNCATION ERROR IS ABOVE THE UPPER BOUND/ DECREASE STEP LEN 

70 NARK«1 

IF (DX-DXMN) 65/65/71 

71 X-X-DX 

C HALVE STEP LENGTH. 

DX*AMAX1(DX2,DXHN) 

C IN CASE MEND HAS BEEN SET/ RESET IT. 

MEND«0 
GO TO 15 

C IF BOTH TRUNCATION ERRORS ARE BELOW THE LOWER BOUND/ INCREASE STFP LEN 

80 IF (DX-DXMX) 81/65/65 

81 X*X-DX 

C DOUBLE STEP LENGTH. 

DXxAMIN1(DX*2./DXHX) 

GO TO 15 
END 


SUBROUTINE P I ( K/ N/ X/ V / XI/ VI) 

C INTERPOLATES BY POLYNOMIAL FITTING. 

C GIVEN V(K)/V(K+l)/*../V(N) AT X ( K )/ X (K+1 )/.../ X ( N > / FITS A SECOND 
C ORDER POLYNOMIAL TC THE THREE POINTS NEAREST XI AND RETURNS A VALUE VI 
C AT XI. FAILS IF ANY TWO OF THE X ARE EOUAL. IF N-K EQUALS ZERO OR ONE 
C THE SUBROUTINE ALSO RETURNS A RESULT. 

C 

DIMENSION X(l)/V(l) 

CHECK IF ONLY ONE DR TWO POINTS IN ARRAY. 

IF (N-K-1) 2/4/6 
C RETURN CONSTANT VALUE IF N*K 
2 V1«V(K) 

60 TO 100 

C INTERPOLATE LINEARLY IF N«K+1 
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4 V1»V(K) + (V(N)-V(K) )/(X(N)-X(K))'MXl-X(K)) 

GO TO 100 

LA6RANGIAN INTERPOLATION, 

FIND NtAREST THREE POINTS. 

6 DO 10 J«K>N 

IF (Xl-X(J)) 15»10»10 
15 I«J 

GO TO 17 
10 CONTINUE 
I«N 

17 IF (I.GT.K+1) GO TO 20 
USE FIRST THREE POINTS. 

I*K-1 
GO TO 25 

20 IF (I.NE.N) GO TO 22 
USE LAST THREE POINTS. 

I*N-3 
GO TO 25 

USE TWO POINTS WHICH STRADDLE XI AND SELECT THE THIRD. 

22 IF (X(I-2)+X(I+l)-2.*Xl) 23#23#2A 

23 I-I-2 

60 TO 25 
2A I«I-3 

INTERPOLATE USING LAGRAN6ES FORMULA. 

25 S«0. 

DO AC J>1»3 
P-1. 

DO 30 L-l»3 
IF (J.EQ.L) 60 TO 30 
P-P*(X1-X(I+L) )/ (X(I+J)-X(I+L) ) 

30 CONTINUE 
AO S«S+V(1+J)*P 
Vl-S 

100 RETURN 
END 


SUBROUTINE SMOOTH(K# L* H,X»Y,DY) 

INPUT Y(I) AS A FUNCTION OF X(I) FOR I-L, L+1# . . . * M. 
SMOOTH IN K PASSES. 

OUTPUT Y(I) SMOOTHED AND THE DERIVATIVES DY(I). 
ROUTINE FAILS IF ANY TWO X EQUAL 
DIMENSION X(1)»Y(1)»DY(1) 

N IS THE NUMBER OF POINTS TO BE SMOOTHED. 

N-M-L+1 

00 NOT SMOOTH IF LESS THAN FOUR POINTS. 

IF (K.EO.O.OR.N.LT.A) GO TO 50 

o 

SMOOTH 

M3-M-3 
DO AO J-l.K 
C SMOOTH SUCCESSIVE SETS OF FOUR POINTS. 

DO AO I-L.M3 
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A«(Y(I + ?»-Y(I) )/(X(I + 2)-X(D) 

C«<Y(I + 3)-Y(I+in/(X(I + 3)-X(l+l) ) 

B«Y(I)-A*X(I) 

D«Y(I+1 J-C*X(I+1) 

C DO NOT SMOOTH IF THE FOUR POINTS ARE COLINEAR. 

IF (A-C) 10#12»10 
12 IF (B-D) 

C XINT IS INTERSECTION OF LINES JOINING ALTERNATE POINTS 
10 XINT»<B-D)/(C-A) 

IF (XINT. GE.Xd+l). AND. XINT. LE.X(I42)) 60 TO 40 
14 Y(I+l)-0.5*(Y(r+l)+A*X(I+l)+B) 
Y(I+2)«0.5*(Y(I+2)+C*X(I+2)+D) 

40 CONTINUE 

CALCULATE DERIVATIVES 

50 IF (N.EQ.l) 60 TO 70 
IF (N.E0.2) 60 TO 80 
C AT FIRST POINT. 

D1*X(L+1)-X(L) 

D2-X(L+2)-X(L> 

51- (Y(L+1)-Y(L))/D1 

52- (Y(L+2)-Y(L))/D2 
DY(L)«(D2*S1-D1*S2) /(D2-D1I 

C AT INTERMEDIATE POINTS. 

Ll-L+1 

M1«M-1 

DO 60 I*L1>N1 
01*X{I)-X(I-1) 

D2«X(I + 1)-X(H 
S1*(YII)-Y(I-1))/D1 
S2=( Y( I+U-Y(n ) /D2 
60 DY(I )=(D1+S2+D2*S1)/(D1+D2) 

C AT LAST POINT. 

D1-X(M)-X(M-1) 

02*X(M)-X(M-2) 

S1«(Y(M)-Y(M-1))/D1 

S2«(Y(M)-Y(M-2))/D2 

DY(M)»(D2*Sl-Dl*S2)/(02-Dl) 

66 RETURN 
C 

C ONLY ONE POINT GIVEN. 

70 0Y(L)*0. 

GO TO 66 

C ONLY TWO POINTS GIVEN. 

80 0Y(L)«(Y(M)-Y(L))/(X(H)-X(L)) 

OY(M)»OY(L) 

60 TO 66 
END 


SUBROUTINE DERIV 

COMMON/XCHGE/ST6N1»P>FF,TORR 

C0MM0N/H0NE/H1>H2 

COMMON/ARRR/KR.NRf XIR(200)f AIR (200 )>DIR (200)f XCR(200) 
COMHON/BBDD/AJ.AAME(200)«AAR(200) 
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COMMON/ ARRME /KME^NME / XIME( 200) ^ AIM E( 200), DIME ( 200) /XCME (200) 
COMMON/ ANS/THT A, DELSTA,GASQT,BHU,GAM21 
COMMON/REIN/TO,TW,5IGFA,GAMMA,REYNO,ASEA 

C0MMCN/C0NST/R,S,G,PI,GW,T0W,B,E91, E92,H,AM£,DME,PT3,REY3f TE,IRY 
COMMON/ INTEG/CF, DEL, DCF, DDL, DX, DXMN, DXMX, X, XEND, MEND 
CALL PI(KME,NME,XCME, AIME,X,AME) 

CALL ?I(KME,NME,XCME,DIME,X,DME) 

CALL PI(KR,NR,XCR,AIR,X,RR) 

CALL PI(KR,NR,XCR,DIR,X,DR) 

SIGNl— !• 

IF(CF.GT.O.) SIGN1»1. 

SIGN»1. 

051*5.1-0.614/ (0.A*ASEA) 

C52*5.1 

G1K*0. 

G2K«0. 

HG1K*0. 

HG2K»0. 

G1*0. 

G2-0. 

HGl-0. 

HG2-0. 

DGlCF-0. 

DG2CF-0. 

DG1DL*0. 

DG2DL-0. 

DGlME-0. 

DG2ME-C. 

HDGlCF-0. 

HDG2CF*0. 

HDG1DL*0. 

HDG2DL-0. 

HDG1ME*0. 

HDG2ME-0. 

AME2»AME*AME 

S1*1.+S*AME2 

C11«GW*S1 

C22*S1-S1+GW 

C33*-S+AME2 

TE*TC/S1 

T0TT«(TC/TE)*+GAH21 

TE198*(TE+198.6)/(T0+198.6) 

REY3=GASQT+TOTT*TE198/BMU 
IFdRY. .LT.l) GO TO 40 
REYN0*PT3+AHE*PEY3/12. 

40 CONTINUE 

ASC*S+AME2*T0W/S1 
A*S0RT( ASQ) 

ASB*2.*ASQ-B 

XY21*2.+S+AHE 

STWE*SQRT(0.5*TW/TE) 

B24=B*B+4.*ASQ 

SQBA*$0RT(B24) 

SARC*ASIN( AS6/SQDA) 

USTT=A/SARC 

USTA*STWE*USTT 

CF3»ABS(CF) 



CF2-SCRT(CF3) 

1FCCF.lt. 0*) SIGN—1. 

SCF*SIGN+CF2 

Rl-SCF+USTA 

Q1*0.5*SQBA/ASQ 

Q2*0.5+B/ASQ 

TETta-TE/TW 

TAT-REYN0+STWE+TETW**1.76 

TATD*TAT^CF2 

TATDC»TATD^DEL 

IFCCF3.lt. TORR) GO TO 83 

R2»2.5+AL0G(TAT0C)+C51 

P»G.5-0.5*R1%R2 

GO TC 8^ 

83 CONTINUE 
P*0.5 

84 CONTINUE 

UU-SQRTCGAMMA*R+TE)+AME*12.+SQRT{G) 

SIG3«S+AME2*SIGMA**0.333333+1. 

TR=T0*SIG3/S1 
GAMl-CGAMMA-1. )*AME 
C DIFFERENTIATE WITH RESPECT TO ME 
DTEM£»-T0*2.+S*AME/(S1*S1> 

0AME«SQPTCT0W*S/S1)*C 1 .-$♦ AME2/S1 ) 

DASQME*2.*A*0AME 

DT17ME«1.76*CTETW)**0.76*DTEME/TW 

DBAME=-0.5*B+DAS0ME/CASQ*ASQ) 

DSQBME«DASQME/ CASQ*5QBA)-0.5*SCBA*DAS0ME/CASQ*ASQ) 

DA2BME=C2.+DAS0ME-C2.*ASQ-B)+2.+DAS0ME/B24)/SQBA 

TIN«1.-CASB*ASB/B24) 

TINSQ*SQRTCTIN) 

DARCME«DA2BME/TINSQ 

DSQTME=-0.5+STWE*DTEME/TE 

DAARME*DAME/SAPC-A*DARCME/(SARC*SARC) 

IFCCF3.LT.T0RR) GO TO 85 

DR1ME»SCF*CDS0TME*USTT+STWE*DAARME) 

DR2ME=2.5+DT17ME/TETW*+1.76+2.5*DSQTME/STWE 

0PME«-0.5*R2*0R1ME-0.5*R1*0R2ME 

DR1CF-0.5+R1/CF3 

DR2CF*1.25/CF3 

DPCF»-0.5+DRlCF*R2-0.5+DR2CF*Rl 
DR2DL-2.5/DEL 
DPDL--C.5+2.5+R1/DEL 
GO TC 86 

85 CONTINUE 
DPME-0. 

DPCF-0. 

DPDL-0. 

86 CONTINUE 

E605=EXPCC10.-C52)/2.5) 

DELlO-OEL/10. 

Y1=E605/TATD 

IF(YI.GT.DELIO) Yl-OELIO 
DY*CDEL-Yl)/50. 

M*1 

N*1 

77 DO 1 I-M#N 
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Y*YX+(51-I)+DY 

YO-Y/OEL 

PYD«PI+YD 

YDDC*YD*+ASEA 

YDD1«1#-YDDC 

YDSQ*SORT(YDDl) 

YDSl-YDSQ+1. 

ALGY*AL0G(YD)+2**C YDSQ-ALOG(YOSl) )/ASEA 
UUBST*1,+2*5*R1«ALGY -P-P+COS ( PYD) 

UUE*Q1*SIN<UUEST*SARC )+Q2 

UUE2-UUE+UUE 

Z«C11+C22*UUE+C33*UUE2 

Zl*l./Z 

Z2«Z1/Z 

B70*CUEST^SARC 

DUSTME«2.5*SCF*ALGY ♦ ( DSQTME*USTT+STHE*DAARME ) -DPME* ( 1 . +CO^ ( P 

lYD)) 

DSINME»C0S(B70)*(SARC*DUSTME+UU£ST*DARCME) 

DUUEME»DSQBME*SIN(B70)+Q1*DSINME+DBAMF 

DZM*(GAMMA-1.)+AME 

DZME*DZM+GW+DZM+(1*-GW)+UU£+DUUEME*C22-DZM*UUE2+2.*C33*UUE*0UU 

lEME 

DUE$TU»01*CDS(B70)*SARC 

DZUUE*C22+2.*C33+UUE 

DUSTCF«2*5*USTA*0,5*ALGY /S CF-D PCF+ ( 1 .+COS ( PYD) ) 

DUUECF»DUESTU*DUSTCF 

DZCF*DZUUE+DUUECF 

RPI-YDSQ/DEL 

DUSTDl«-2.5*USTA*SCF*RPI-DPDL*(l,+COS( PYD) >-P*SIN{PYD)+PYn/DEL 

DUUEDL*DUESTU+OUSTDL 

DZDL»DZUU£*DUUEDL 

G1KP*UUE 

G2KP«UUE2 

GIP^UUE+Zl 

G2P-UUE2+Z1 

DG1CFP»Z1*DUUECF-Z2^UUE*DZCF 

DG2CFP*^2. + UUE*CUUECF*Z1-UUE2*DZCF + Z2 

DG1DLP«Z1*DUUEDL--Z2 + UUE*DZDL 

D62DLP*2.*UUE*DUUEDL+Z1-UUE2+DZDL^Z2 

DG1MEP=Z1*DUUEME-Z2*UUE*DZME 

DG2H£P*2.+UUE*DUUEME*Z1-UUE2*DZME*Z2 

G1K«G1K+G1KP 

G2K«62K+G2KP 

Gl-Gl+GlP 

G2«G2+62P 

DG1CF=DG1CF+DG1CFP 
DG2CF*DG2CF+DG2CFP 
DG1DL»DG1DL+DG1DLP 
DG2DL=DG2DL+DG2DLP 
DG1ME=DG1ME+DG1MEP 
D62ME«DG2ME+DG2MEP 
1 CONTINUE 

HG1K=HG1K+G1KP 

HG2K*HG2K+62KP 

HG1«HG1+G1P 

HG2-HG2+G2P 

HDG1CF*H0G1CF+DG1CFP 
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hDG2CF»HDG2CF+DG2CFP 

HDG1DL*HDG10L+DG1DLP 

HDG2DL*HDG2DL+062DLP 

hDGl«E*HOGiriE + DGlMEP 

H0G2ME*HDG2ME+0G2MEP 

IF(M.NE.l) 60 TO 11 

M«2 

N*51 

60 TO 77 
11 CONTINUE 

G1I*(G1“0.5+HG1)*DY+0.5*Y1*61P 

G2I»(G2-0.5+HG2)*DY+0.5*Yl*G2P 

G1KI=(G1K-0.5*HG1K)*DY+0.5*Y1*G1KP 

G2KI*(G2K-0.5*HG2K)*DY+0*5*Y1+G2KP 

THTA*G1I-62I 

DELSTA«DEL-G1I 

DG1CFI»(DG1CF->0.5*HDG1CF)*DY+0.5*Y1^DG1CFP 

DG2CFI*(D62CF-0.5*H0G2CF)*DY+0.5+Yl*0G2CFP 

DGlDLI*CDGlDL-C.5*HDGiDL)*DY+C.5*Yl+DGlDLP 

DG2DLI-(DG2DL-0.5*HDG2DL)*DY+0.5+Yl*0G2DLP 

DGlHEI*(DGlME-0.5^HDGlPE)*0Y+0.5+Yl*DGlMEP 

DG2MEI={DG2ME-0*5*H0G2ME)*DY+0.5*Y1*DG2MEP 

DTHCF*DG1CFI-DG2CFI 

DTHDL»DG1DLI-062DLI 

DTHME*DG1«EI-DG2MEI 

DDSCF*-D61CFI 

DDSM£«-DG1MEI 

DDSDL — DGIDLI 

H*0ELSTA/THTA 

A8*( (H+1.-TP/TE)*TETW-1. )/l,12 
IF(A8.LT.O.) WRITE(6#700) A8#H 

700 F0RMAT(5X>3HAe«,F10-6>12X>2HH«fF10.6) 

IF(A8.LT.0.) WPITE(6#701) GlI# G2I# DEL> DELSTA^ CF^ DX 

701 FORMAT< 5X,4HG1I«,F10.6>5X, AHG2I*>F10.6#5X, AHDEL=^F10.6#'?X# 
17H0ELSTA«,F10.6#3X#3HCF->Fi0.6>3X,3hDXe,F10.6) 

AA«A8**E91+2. 

A1-*E91*A8+*E92 /1. 12 

DTRDX*DME+GAMl+T0*(SIGMA++.333333-TR/T0)/$l 
FRAC*(4.-2. + AA)*(A.--2.^AA)/(2.-8. + AA+2.*AA*AA) 

Hl«( 1,-AA+AA)/ (4.-2.+AA) 

H2»DEL/THTA-H 

H2K*G1KI/(G1KI-G2KI) 

IF(CF.GT.O.) GO TO 50 


FF«(H2+A,2)*,0067 
GO TO 51 

50 CONTINUE 

FF*0 •0306* (H2K-3.0) ♦♦-0.653 

51 CONTINUE 

DUDX1*12.+SQRT(G+6AMMA+R) 

TES2«T0^AME2^(GAMMA-1. )/(2.+SQRT(TE)^$l+Sl ) 
0UDX»DUDX1^(SQRT(TE)-TES2)^DME 
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DUDXU«THTA1>DUDX/UU 

C 

DH20X«(FF-H2*CF/2. )/THTA+H2*(H+L*»*OO0*/aO: 

DTHDX*CF/2.-(H+2.-AME2)*0ti0XU 

IF(AJ.NE.O.>DTHDX*DTHOX-THTA*DR/l» 

0TEDX«-TC*2.*S*AME*0ME/{Sl*il) 

T15»(DTRDX-DTEDX)*THTA/TE 

DDSD1*DDSDL-1. 

DETERM*DTHCF*0DSD1-DTHDL*0DSCF 
THME«DTHDX-DME*DTHME 
DDSDC»-THTA*DH20X-H2*0THDX 
DSME«D0SDC-D»1E*0DSME 
0CF«(0DSD1*THME-DTHDL*DSME)/DETERH 
DOL-(DTHCF*DSME-DDSCF*THME)/OETERM 
DOSDX-DDL+DDSDC 
30 RETURN 
END 


TABLE 5-B: INPUT TO PROGRAM BL6RN 


2.82 10 + 

13 TEETER 





536.6 

536.6 

73. 

0.72 

.00169 

0.126 

65 

0.0001 

0.01 

1. 

1. 



1 

1 . 







4. 






<r514 

2.4427 

2.4187 

2.403 

2.3896 

2.3802 

2.3648 

65 

2.708 

4. 

2.814 

2.886 

2.9542 

3.0018 

3.08 


2.65 


2.706 2.81A 2.886 2.95A2 3.0018 3.08 


100 



TABLE 6 


INPUT TO METHO D OF CHA RACT ERIST ICS PROGRAM 


10^13 

FROM 2.4 




2.65 





0 1 

k7 1 

0 0 

0 


2.62 

1.4 

24. 

25.31 

4.1 

0.0 COCCI 

0.001 

0.00001 

0.001 


3.00 

0.001 

0.001 

0.0 

5. 


3. 




1.942 

3.235 

3.25 

5. 


0.342375 





13. 

13. 

-10. 

-10. 


2.02 

3.02 

5.0 



1.02 





0. 

0. 

0. 



6.335 

0.999 

2.584 



2 

3 2 

1 2 : 

26 


2.40 

1.02 

0.0 

2.53637 

0.99635 

2.40 

1.01 

0.0 

2.53637 

0.99835 

2.40 

1.00 

0.0 

2.53637 

G.99835 

2.40 

0.99 

0.6 

2.53637 

C.99835 

2.4G 

0.9797 

1.2 

2.5369 

G.99835 

2.40 

0.971627 

1.8 

2.538211 

Q. 99835 

2.40 

0.9617 

2.4 

2.66272 

0.999 

2.40 

0.9516 

3.0 

2.6588 

0.999 

2.4C 

0.9416 

3.631 

2.655 

0.999 

2.40 

0.9315 

3.7305 

2.6515 

0.999 

2.40 

0.921477 

3.8277 

2.6482 

0.999 

2.40 

0.91144 

3.923 

2.645 

0.999 

2.4« 

0.9014 

4.0164 

2.6418 

0.999 

2.40 

0.861519 

4.19945 

2.6359 

0.999 

2.40 

0.86174 

4.3785 

2.63035 

0.999 

2.40 

0.84213 

4.555 

2.62511 

0.999 

2.40 

0.8227 

4.72992 

2.620 

0.999 

2.4C 

0.80 348 

4.904 

2.6154 

0.999 

2.40 

0.78445 

5.078366 

2. 6109 

0.999 

2.40 

0.765465 

5.25314 

2.60662 

0.999 

2.40 

0.7471 

5.4289 

2.6025 

0.999 

2.40 

0.728777 

5.6062 

2.598 

0.999 

2.4C 

0.7J.37 

5.7853 

2.5948 

0.999 

2.40 

0.69297 

5.96648 

2.5912 

0.999 

2.40 

0.67553 

6.15 

2.58775 

0.999 

2.40 

0.6584 

6.3367 

2.58447 

0.999 

1 

3 3 

1 1 

9 


2.40 

0.44812 

13.0 

2.42914 

C.998 

2.40 

0.466C1 

12. 33 

2.4298 

0.998 

2.40 

0.50812 

11.2254 

2.43212 

0.998 

2.40 

0.51993 

10.9846 

2.43446 

0.998 

2.40 

0.581397 

9.6785 

2.447921 

0. 998 

2.40 

0.60356 

9.5647 

2.451 

0.998 

2.40 

0.616255 

9.44656 

2.4505 

0.998 

2.4 

0.624544 

9.375146 

2.4500 

0.998 

2. 4012 

0.63984 

9.24714 

2.4495 

C.998 
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ooonooonnoo oooooo 


TABLE 7-A; PROGRAM MFLX 


PROGRAM HFLXIINPUT,OUTPUT*TAPE 5=INPUTt TAPE 6=0UTPUTI 


PROGRAM FOR COMPUTATION OF MASS AND MOMENTUM FLUX 
INPUT FORMAT 7F10.6 EXCEPT CARO 1 
CARO (SI COLUMNS 

1 TITLEf COLUMNS 1-72 HOLLERITH 

2 1-10 TO TOTAL TEMPERATURE DEGREE R 

11-2C CR GAS C0NSTANT(=53.3I 

21-30 AN NO. OF POINTS INPUT 
31-40 DY AVERAGE INCREMENT (INCH) 

Ll-5 0 RC RADIUS OF LOCAL CONE SURFACE CINCH) 
51-60 AIN =0. NO MORE JOB AFTER THIS INPUT 
=1. MORE JOB AFTER THIS INPUT 

3 1-73 R RADIUS AN VALUES CINCH) 

4 1-70 EM MACH NUMBER AN VALUES 

5 1-70 P STATIC PRESSURE AN VALUES CPSIAI 

DIMENSION EMC2Q0) »R( 200 ) , ROU ( 20 0 ) ,ROUUC200) 
DIMENSION PC200) 

UTMENSION TITLEC12) 

1 REAO(5*1000) (TITLEC J)» J=l»12) 

WRITE (6,20 00) ( TITLE (Jl , 1* 12 ) 

REAO(5,1QO) TO,CR,AN,OY,RC,AIN 

WRITE(6,400I TO,RC 

N=AN 

READ (5f 100) (R(I) ,I=1»N) 

READ (5,1 0 01 (£M(I) ,I=1,N) 

READ (5, 13 0) (Pdl ,1=1, NJ 
WRITE(6,300) 

SUM=0. 

BU=0. 

BUU=0. 

SUN=0. 

DO 2 1=1, N 
AI=I ^ 

A=SQRTCTOI 

TOT=l .♦0.2»EM(I)»EM(I) 

A3=SQRT(TCT) 

R0U(I)=4 9.-»144.»P(I)^EM(I)»AB/(CR^A| 

U=49.^EM( II^A/AB 
RCUU(I)=RCU(I)^U 
IF(I.EQ.l) GO TO 20 o 

RIU=0 • 5*i ROU ( I-ll +ROU Cl ) ) 

RIUU=0.5» (ROUU(I-l)+ROUU(I)l 

AR£=5.1416» (R(I)*R(I)-R (1-1)*RC I-l))/144. 

BU=RIU»AR£ 

BUU=RIUU*AR£ 
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20 SUM=SUM*8U 
SUN=SUN4-BUU 

URITEI6«Z00) Rdl «P ( li t SUH, SUN 

2 CONTINUE 
1000 F0RMAU12A6) 

2000 FCRNAT(lHl,lGXtl2A6//l 
ICO FORMAT <7F10.6> 

200 FCRMATC 10X«Fi0 .6 *SX *F10. 6* 5X *F 10 .6* 5X « FIO . 6,5X t F15 .6) 

300 FCRHATJ/15X,lHR»15XtlHM,15X,lHP,lCX,9HMASS FLUX»10X,8HMOH FLUX//) 
400 FORHAT(/20X,3HTO=tF6.2t9H DEGREE R tSXtAH RC= vFS.4t4H IN./) 
IFfAIN.EQ.l.) 60 TO 1 
END 


TABLE 7-B: INPUT TO PROGRAM MFLX 


M=2.82 

X=2.85 

MASS FLUX 

CHECK 




536.6 

53.3 

70. 

0.005 

0.55145 

1. 


0.551i^5 

0.555 

0.560 

0.565 

0.570 

0.575 

0.580 

0.565 

0.590 

0.595 

0.600 

0.60 5 

0.610 

0.615 

3.620 

0.625 

0.630 

0.635 

0.640 

0 .645 

0*650 

0.655 

0.660 

0.665 

0.670 

0.675 

0.680 

0.685 

0.690 

0.695 

0.7CQ 

0.705 

0.710 

0.715 

0.720 

0.725 

0.730 

0.735 

0.740 

0.745 

0.750 

0.755 

0*760 

0.765 

0.770 

0.775 

0.780 

0.785 

0.790 

0.795 

0.800 

0.805 

0.810 

C.815 

0.820 

0.825 

0.830 

0.835 

0.840 

0.845 

C.850 

0.855 

0.860 

0.865 

0.870 

0.875 

0.880 

0.885 

0.690 

0.895 

2 .if 34 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2.434 

2. '43 4 

2.434 

2.434 

2.43 4 

2.434 

2.434 

2.440 

2.440 

2.440 

2.440 

2.440 

2.440 

2.440 

2.440 

2.440 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.450 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.460 

2.400 

2.335 

2.335 

2.335 

2.335 

2.335 

2.32 5 

2.325 

2.320 

2.310 

2.310 

2.311 

2.340 

2.347 

2.345 

2.345 

1.810 

1.813 

1.610 

1.810 

1.810 

1.810 

1.810 

1.8C9 

1.810 

1.610 

1.805 

1.805 

1.805 

1.805 

1.805 

1.8Q5 

1.800 

1.800 

1. BOO 

1.800 

1.800 

1.795 

1.795 

1.795 

1.795 

1.790 

1.790 

1.790 

1.785 

1.785 

1.780 

1.780 

1.770 

1.775 

1.775 
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1.770 

1.770 

1.765 

1.760 

1.755 

1.754 

1.750 

1.750 

1.740 

1.734 

1.723 

1.720 

1.720 

1.720 

1.725 

1.725 

1.735 

1.740 

1.630 

2.110 

2.120 

2.150 

2.160 

2.170 

2.170 

2.170 

2.170 

2.170 

2.170 

2.170 

2.150 

2.130 

2.125 

2.090 

2.030 


X=2. 

536.6 

968 MASS 

53.3 

CHECK 

69. 

9.005 

0.58331 

1. 


0.56331 

0.585 

0.590 

0.595 

C.60C 

0.6C5 

0.610 

0.615 

0.620 

0.625 

0.630 

0.635 

0.640 

0.645 

0.650 

0.655 

3.660 

0.665 

0.670 

0.675 

0.680 

9.685 

0.690 

0.695 

0.709 

C.705 

0.710 

8.715 

0.720 

0.725 

0.730 

0.735 

0.740 

0.745 

3.750 

0.755 

0.76Q 

0.765 

0.770 

0.775 

Q.78Q 

0.785 

0.790 

0.795 

0.800 

0.805 

0.810 

0.815 

0.620 

0.825 

0.830 

0.835 

0.340 

0.845 

0.850 

3.855 

0.860 

0.865 

0.670 

0.875 

0.880 

0.865 

0.890 

0.895 

0.908 

0.905 

0.910 

0.915 

0.920 


2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2.444 

2. 307 

2.337 

2. 307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.307 

2.3C7 

2.307 

2.307 

2.307 


1.78*1 

1.784 

1.784 

1.764 

1.784 

1.784 

1.784 

1.7 84 

1.784 

1.784 

1.784 

1.784 

1.784 

1.784 

1.784 

1.784 

1.764 

1.784 

1.784 

1.784 

1.764 

1.7 64 

1.784 

1.784 

1.784 

1.784 

1.784 

1.784 

1.785 

1.785 

1.785 

1.785 

1.795 

1.795 

1.795 

1.795 

1.796 

1.796 

1.798 

2.220 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.215 

2.220 

2.220 

2.220 

2.22C 

2.229 

2.220 

2.220 

2.220 

2*220 

2.220 

2.220 

2.220 

2.220 

2.220 

2.22 0 

2.220 

2.220 


X=3.205 

MASS 

FLUX CHECK 





536.6 

53;3 

63. 

0.005 

0.6334 



0.633417 

0.635 

0.643 

0.645 

0.650 

0.655 

0.660 

0.665 

0.678 

0.675 

0.680 

0.685 

0.690 

0.695 

0.7C0 

0 .705 

0.710 

0.715 

0.720 

0.725 

0.730 

0.735 

0.740 

0.745 

0.759 

0.755 

0.760 

0.765 

0.779 

0.775 

0.780 

0.785 

0.790 

0.795 

0.800 

3.8C5 

0.610 

0.815 

0.820 

0.825 

0.830 

0.835 

0.840 

0.845 

0.650 

0.855 

0.860 

0.8 65 

0.870 

0.875 

0.880 

0.885 

0.890 

0.695 

0.900 

G.905 

0.910 

0.915 

0.920 

0.925 

0.930 

0.935 

0.940 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.445 

2.350 

2.295 

2.260 

2.265 

2.250 

2.240 

2.245 

2.250 


104 



I 


2*2 55 

2.260 

2.270 

2.283 

2.300 

2.310 

2.315 

2.315 

2.310 

2.305 

2.309 

2.309 

2.360 

2.380 

2.375 

2.370 

2.265 

2.262 

2.260 

2.260 

2.260 

2.256 

2.255 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

2.252 

1.775 

1.775 

1.775 

1.775 

1.77 5 

1.775 

1.775 

1.775 

1.775 

1.775 

1.775 

1.775 

1.775 

2.000 

2. 2^*7 

2.250 

2.260 

2.275 

2.298 

2.290 

2.280 

2.270 

2.269 

2.250 

2.240 

2.230 

2.220 

2.220 

2.182 

2.192 

2.21D 

2.230 

2.25G 

2.270 

2.290 

2.310 

2.335 

2.350 

2.365 

2.370 

2.375 

2.380 

2.384 

2.399 

2.395 

2.400 

2.401 

2.401 

2.401 

2.401 

2.401 

2.401 

2.401 

2.401 

2.401 

2.410 

2.410 

2.410 

2.415 

2.415 

2.415 

2.415 

2.415 


TABLE 8: INPUT TO PROGRAM ANAL 


SECOND SHOCK M=2.82 10+13 

SLOT BLEEO=0.028 


2.411 

0. 037747 

77599.4 

2. 

1. 

-0.028 

0. 

2. 

0.4 

5.1 

4.5 

33. 




0.1225 

1.0 

0.5 



13. 





1.00 

1. 00 

2.3447 

1.156 

4.52 

1.03627 

1.005586 

2.322 

1.1987 

4.02 

1.072559 

1. 00 9497 

2.32 

1.2014 

3.998 

1.1D013 

1.0134C7 

2.2993 

1.24084 

3.954 

1.12697 

1.01676 

2.293 

1.25415 

3.90 

1.18140 

1.0229 

2.269 

1.3007 

4.00 

1.21767 

1.027S33 

2.259 

1.32075 

4.010 

1.25758 

1. 0305 

2.248 

1.3457 

4. 018 

1.27935 

1.03575 

2.244 

1.35682 

4.02 

1.29749 

1. 038 

2.24 

1.3667 

4.04 

1.36279 

1. 0447 

2.234 

1.40843 

4.08 

1.39907 

1. 05027 

2.230 

1.44133 

4.12 

1.43535 

1. 05586 

2.227 

1.48133 

4.16 


NASA- Langley, 1976 
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